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OBJECTIVE 


Translate well-known differential equation solutions into a working program to com- 
pute propagation in underwater acoustic ducts. Document the program methods, to assist 
users of this and similar programs. 


RESULTS 


1. An effective program for computing propagation loss in a layered ocean by normal 
modes has been developed. Complete documentation for the program is contained herein. 


2. Sediment layers are modeled as fluids in which densities, sound speeds, and absorp- 
tion can be specified. This permits a complete wave solution for bottom reflected sound 
energy. 

3. A continued fraction technique for evaluating asymptotic series is shown to give 
superior results in evaluating the auxiliary functions required in this program, the modified 
Hankel functions of order 1/3. 

4. A mode follower program given here is useful in tracing eigenvalues. Such traces 
are needed to understand the eigenvalue structure. 


RECOMMENDATIONS 


1. Improve the mode locating ability of this normal-mode program to make it self- 
contained. It currently requires user interaction to locate eigenvalues. 

2. Investigate methods to incorporate the effect of rough boundaries into this 
program. 
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INTRODUCTION 


This report describes a normal-mode program that has been used successfully for 12 
years to compute sound propagation in idealized underwater acoustic ducts. The theory and 
considerations used in developing the program are discussed here, and a copy of the FOR- 
TRAN statements are included as appendix A. Appendix B consists of sample inputs and 
outputs to assist users in gaining familiarity with the program. It is hoped that this report 
contains sufficient information to allow a user to run the program and to modify it as desired. 


This program follows the methods developed by Furry and Freehoffer (ref 1) to com- 
pute electromagnetic propagation in the 1940s. Marsh adapted these methods to underwater 
sound in his doctoral thesis (ref 2). Using this material, Pedersen, at NOSC in the late 1950s, 
adapted the method to digital computers and developed the programs to compute the auxil- 
iary functions. This original program used two layers to define the sound-speed profile (ref 3). 
This program was expanded to three layers by DF Gordon and RF Hosmer and finally to the 
multiple-layer program reported here. In this program the only constraints on the number of 
layers are computer space and running time. The program is normally configured to permit 
up to 12 layers. 


The earlier programs were used to study sound propagation in ocean surface ducts. 
Programs that permit more layers have proven useful also for studying propagation in the 
deep ocean, although the number of modes required generally limits computations to fre- 
quencies below 300 Hz. The multiple-layer program has also proven useful in modeling sedi- 
ment layers and thus in computing shallow-water propagation. 


The principal limitation in the application of this program to real-world situations is 
the requirement of ideal conditions: boundaries must be smooth and horizontal, and no 
variation of boundary conditions with range is permitted. Despite this limitation, the pro- 
gram has proven useful in predicting and explaining acoustic propagation and has applications 
in a number of related areas. These include checking other types of wave-theory models or 
corrections such as caustic corrections; determining group velocities, dispersion curves, and 
reflection coefficients; and determining acoustic coupling between ducts. 


The following paragraphs describe the specific topics covered by the sections in this 
report. In GENERAL SOLUTION are the equations required to solve the wave equation 
with the boundary conditions used here. DETERMINANT is part of the basic solution but 
is concerned with the particular numerical method used in this program to evaluate the con- 
ditions imposed by the boundaries. Other approaches could be used instead. A later section, 
NUMERICAL BREAKDOWN, is also part of the basic solution, but deals with special numer- 
ical problems that have arisen but are not apparent from the basic equations. 


1. The Bilinear Modified-Index Profile, by WH Furry, in Propagation of Short Radio Waves, DE Kerr, ed; 
MIT Rad Lab series, vol 13, p 140-168, McGraw-Hill, New York, 1951. 

2. Navy Underwater Sound Laboratory Report 111, Theory of the Anomalous Propagation of Acoustic 
Waves in the Ocean, by HW Marsh, 1950. 

3. Normal-Mode Theory Applied to Short-Range Propagation in an Underwater Acoustic Surface Duct, by 
MA Pedersen and DF Gordon; J Acoust Soc Am, vol 37, p 105-118, January 1965. 


FINDING EIGENVALUES deals with the philosophy of eigenvalue location employed 
by this program, which essentially leaves this function to the user, the program only serving as 
a tool. It shows how the program is used to make computations. 


Several “‘automatic”’ mode finding versions of this program have been developed to 
the point of accommodating certain classes of profiles. However, they need further develop- 
ment and have not yet been reported. 


SOUND SPEED PROFILE indicates the required equations for curve fitting and the 
various ways the sound speed can be read in on cards. A continuous water profile can be 
entered quite simply, but sediment layers with sound speed discontinuities and absorption 
gradients can become complicated. 


REFLECTION COEFFICIENTS AND OTHER AUXILIARY OUTPUTS describes a 
short subroutine that computes reflection coefficients for any mode at a given profile inter- 
face. Intermode interference lengths and mode damping coefficients are also discussed. 


COMPUTATION OF THE MODIFIED HANKEL FUNCTIONS gives the analysis 
necessary for computing these functions. The use of continued fractions to evaluate an 
asymptotic series is discussed. To facilitate running the program on computers of different 
word length, this section provides the information required to optimize the functions for 
the different word lengths. 


MODE FOLLOWER PROGRAM describes a separate but related program for investigating 
the eigenvalues themselves rather than using them to compute propagation losses. 


GENERAL SOLUTION 


The derivation of the normal-mode solution has been discussed from various points 
of view (eg ref 1, 4,5). Only an outline is given here. In general, the time-independent wave 
equation is written in polar coordinates and the azimuthal coordinate is dropped under the 
assumption that the field is independent of azimuthal direction. Thus. 


(1/r) (8/8r) [r(8W/dr)] + (02/822) + (w2/c2) Y =0, (1) 


where w is the velocity potential, c the sound speed, and the independent variables are depth, 
z, and range, r. 


Equation (1) is then separated into range- and depth-dependent parts with a separation 
constant A’. The separation is possible when the sound speed is a function of depth only. 
After accounting for the source discontinuity and the outgoing radiation condition, integrating 
over all real values of the separation constant, and normalizing, one can find the solution for 
a field point in terms of propagation loss H as follows: 


4. Naval Air Development Center Report NADC-72002-AE, Normal Mode Solutions and Computer Programs 
for Underwater Sound Propagation, by CL Bartberger and LL Ackler, 4 April 1973. 

5. A Normal Mode Theory of an Underwater Acoustic Duct by Means of Green’s Functions, by RL 
Deavenport; Radio Sci, vol 1, p 709-724, 1966. 


N D 
D 
H = -10 log P Ppt », Ho Ant) Un) U,, Zo) + Mat, (2) 
n=1 


where r is the range, zg is the source depth, z is the receiver depth, Ho is the Hankel function 
of order zero, second type, Ap is the nth eigenvalue, Up, is the depth function for mode n, 
and p, and pp are the densities at source and receiver. The sum is over the number of modes, 
N, making a significant contribution. The final term contains the volume attenuation coeffi- 
cient, aA. From Thorp (ref 6), vA in dB/m is computed by the relationship 


0.9144 a, = 0.0001 F2/(1 + F*) + 0.04 F7/(4100 + F?), (3) 


where F is the frequency in kHz. Improved equations or those for specific ocean areas can 
be easily substituted. The depth function, Up, is a solution to the depth-dependent part of 
the separated wave equation 


d2U/dz2 + [w2/c2(z) - 2] U=0, (4) 
where 
w = 2nf 


and f is the frequency, in Hz. 


A closed-form solution to eq (4) can be obtained when the reciprocal sound speed 
squared or squared index of refraction is a linear function of depth. That form is used in 
this program, and sound speed in each layer is expressed as follows: 


[c/e(z)]2= 1-24; @-z)c;, (5) 


where cj, z;, and yj; are the sound speed, depth, and sound-speed gradient, respectively, at the 
top of layer i. Up to 12 such layers are permitted by the program, for modeling the sound- 
speed profile. 


With this expression for sound speed, solutions to eq (4) can be expressed in terms 
of solutions to Stokes’ equation 


h’ +zh=0. (6) 
Only a simple change in independent variable is required from z to §, where 
3 2 
GO) = E (=e) w/c: - 2 a? (7) 
and 
Bile WO) ys) 
ae -2y; w Ic. : (8) 


6. Analytic Description of the Low-Frequency Attenuation Coefficient, by WH Thorp; J Acoust Soc Am, 
vol 42, p 270, 1967. 


The solutions to Stokes’ equation that are used are the modified Hankel functions of order 
1/3, hy(¢) and h9(f). The depth function is a linear combination of these two independent 
solutions: 


Fy i(2)= Ap ty Gp) + By jhoGy), (9) 


where F,, is the unnormalized form of Un. The coefficients apy j and By j for mode n in layer 
i are determined to satisfy boundary conditions. which will be listed below. Values of A, for 
which the boundary conditions can be satisfied are the eigenvalues. 


The first boundary condition is the radiation condition. It is satisfied by using a nega- 
tive sound-speed gradient in the deepest layer, which extends to infinite depth, and by letting 
the depth function there be proportional to hj only. That is, 


F(Z) = Bho Gy). (10) 
At the surface the depth function is zero: 

ei CO) =O), (11) 
and at !ayer interfaces, o U and its depth derivative are continuous: 

Pin i) = Pit Pp it 1 @): (12) 

dF, j(2)/dz = dF, i+ (2)/dz. (13) 
Here pj is the density in layer i, and the excess acoustic pressure, p, is given by 


p=pu. 


If U is assumed to be the vertical component of the velocity potential, eq (12) and (13) are 
equivalent to requiring that the pressure and the vertical component of particle velocity be 
continuous across the layer interface. 


Applying these boundary conditions to a sound-speed profile consisting of M layers 
results in 2M - | linear equations in hy and hy. They are homogeneous in that the constant 
is zero in each equation. There are M-1 coefficients Aj to be determined and M coefficients 
B;. These coefficients can therefore be determined within a constant of proportionality D, 
provided the system of equations is linearly dependent. That is, the 2M - | square matrix of 
coefficients of A; and B; must be of rank 2M - 2 or less. Its determinant will then be zero. 
This is the eigenvalue condition. Values of \ must be found which make the determinant 
zero. This determinant, G, is discussed in more detail in a later section. 


Zeroes of the determinant, G, are found by using the secant method. The variable in 
this iterative method can as well be some function of \ as Ad itself, and we use the following 
complex phase velocity (v): 


he = w/v 


ne (14) 
To find a v that is a root of G requires an initial guess, vj, where the subscript | refers to the 
step in the iteration and a small increment, 6;. Each succeeding estimate is given by the 
relationship 


Vj+1 = Vj oF 5; 
where 

8) = -(j - V4) GIG; - Gj_}). (15) 
The details of this iterative process are given in a later section. 


When an eigenvalue Wea is found, the coefficients are then evaluated. One coefficient 
can be given an arbitrary value, so Aj is set to ojh9[51(0)]. From eq (11), By is then 
-p hy, [¢1(0)]. Pairs of equations (eg (12) and (13)) for each successive interface can then 
be used to evaluate the next Aj; and B; as discussed later. 


Finally the normalizing factor, Dae for mode n is obtained by the relationship 


Daz J oF2@az. (16) 
0 


This equation follows from the orthogonality of the depth functions. It is not the pressure, 
however, which is proportional to pU, but p?2U that is orthogonal (ref 7). Therefore, D, 
must be determined such that the integral of pU2 is 1. 


From Stokes’ equation (eq (6)) and eq (7-9), the integral of F2 takes the form 


Zi+] 3 Zi+] 
i) F2($)dz = [ne Fe + rei] . (17) 
2 


Zi 
Therefore 


n-1 
3 De 
1D = =; W7/ay + » Oj LS j(Z54 1 )/aj - P4554 1 (Zi 1)/j4] Pj41)1 Fe (Zi41) 
i=] 


3) 3) ' 
+l - Piet!) F741) ’ (18) 


where eq (12) and (13) have been used to combine terms at each interface. The derivative of 
F takes the form 


Fi(Zi4 1) = 3j { Ain (54241) % Bins (54@+1)1| : (19) 


The Wronskian, W, is an imaginary constant (see eq (85)) and is the contribution of eq (17) at 
the surface: 


W = -1.45749544104i. 


7. Some Effects of Velocity Structure on Low-Frequency Propagation in Shallow Water, by AO Williams; 
J Acoust Soc Am, vol 32, p 363-365, March 1960. 


The depth functions are normalized by the relationship 
U, (Zo) Up = Fy(Zo) F,p@/Dp- (20) 


The functions F and F’ used in computing Dy are conveniently assembled from the 
elements of the determinant and the coefficients Aj and B;. This requires care in developing 
the computer code, because F is always multiplied by p and F’ has the term a; in it. The sur- 
face differs from the other layers in that Fy is zero there and Bie by eq (19), isayW. How- 
ever, because p; appears as a factor in the coefficients of Fj, the actual value of Fi at the 
surface in the computation is p] ajW. This factor of p1 together with the ar needed for 
orthogonality, when squared, gives the p; of eq (18). 


DETERMINANT 


Normal modes are determined by finding the eigenvalues of a characteristic equation 
which, in turn, is obtained by setting a determinant to zero. The determinant is obtained 
from the coefficient matrix of a set of linear, homogeneous equations expressing the bound- 
ary conditions as given by eq (10) — (13). Since the method of handling this determinant is 
a central feature of this normal-mode program, it is given in detail here. 


The first line of the matrix is taken from eq (11) as 
By ey ho lS, (0)) + Ay py hy [$y (O)]=0. (21) 


At each profile interface, i, where i numbers the interfaces below the surface from | to 
N-1, the two boundary conditions given by eq (12) and (13) are 


Bj jh (5; (2 j41)] + Aj oj by (8; (Z i44)] - Bi+y Oj4) bo USi4y Zi41)] 
alan Onsey Wn Mee Caen) = O (22) 
and 
Bj aj hy 15; (Zig g)] + Aj aj Dy (65(2141)] - Bigy aigy Do (Sj41 Zip] 


~ Ais] any by in] Zep] = 0 - (23) 


The coefficients of Aj in the first equation and Bj+] in the second will be the diago- 
nal elements of the matrix. The nonzero elements of the matrix will therefore be no more 
than two places from the diagonal. The matrix can be stored in the computer in an array of 
size (2M-1) X 4, where M is the maximum number of layers in the sound-speed profile. In 
the final layer, AN hy is omitted, as in eq (10). In the program, the real and imaginary parts 
are stored in separate arrays. 


The sparseness of the matrix permits efficient evaluation by a triangularization 
process of row reduction. For each pair of rows representing a pair of equations given by 
eq (22) and (23), the first element from the first equation and the first two from the second 
equation must be set to zero by subtracting the proper multiple of preceding rows. The 
determinant is then the product of the diagonal elements of the triangularized matrix. The 
value of the determinant, G, is used in eq (15) to find the roots by iteration. 


Note that a value of v that makes this determinant zero, or near zero, ordinarily is 
zero because only one diagonal element is very small. For trapped modes this element is at 
the row representing the first interface below the mode, ie the interface just below the layer 
of positive gradient in which the sound speed is equal to the mode phase velocity. For 
unstrapped modes it is usually the final diagonal element that is small. Thus the layers in 
which the sound speed is greater than the phase velocity of a mode do not greatly affect the 
eigenvalue. Eigenvalues are determined mainly by those parts of the sound-speed profile 
that are less than the phase velocity. 


When an eigenvalue is found, the coefficients Aj and Bj must next be evaluated. As 
mentioned earlier, one coefficient can be arbitrarily chosen. This is done, and eq (21) is satis- 
fied by letting 


Ay = Py holF, (0)] 
and 


By, =-e 1h, [h,(0)]. (24) 


The factor p1 is used simply because the number containing it is easily available in the pro- 
gram. It is divided out by the normalizing factor, D. Eq (22) and (23) can then be used to 
evaluate the remaining coefficients, but the triangularized form of the matrix yields the coef- 
ficients with less computation. If Sij is the element in the ith row and jth column of the tri- 
angularized matrix, then by Cramer’s rule, 


Bi = Aj] 82)-2, 2-2 82i-1,2i/Fi 
and 
Aj = -Aj_1 82i-2,2i-2 82i-1, 2i-1/ Ei 
where 
Ei = £2i-2, 2i-1 82i-1,2i ~ 82i-2, 21 82i-1, 2i-1- (25) 


A simpler form is used for By in the final layer since there is no Ay there. 


In certain situations numerical problems can arise in evaluating the determinant. 
These require some extra tests in the subroutine that makes the evaluation. The extra tests 
will be discussed in the section, NUMERICAL BREAKDOWN. A more routine problem is 
the loss of accuracy that can arise in subtractions in the row reduction of the matrix. This 
loss results in less sharpness of convergence to aroot. The size of the determinant, G, can be 
14 orders of magnitude less at a root than at the general background near the root. This 
variation occurs because the modified Hankel functions can be computed to about 14-place 
accuracy in a computer with 18 decimal places available. Modes usually converge to 10 or 
12 places; thus a few places are lost in evaluating the determinant. In some profiles, usually 
those with multiple ducts or those in which propagation through bottom sediments plays a 
large part, the convergence can be much poorer. Modes need to converge to about 4 places 
to be reliable for computing losses, and convergence occasionally fails to meet this require- 
ment. The only current cure for this loss in accuracy is to go to higher-precision arithmetic 
or to compute the modified Hankel functions to greater accuracy. For instance, a standard 


matrix triangularization routine that uses full row and column pivoting has been tried with 
no resultant increase in accuracy. 


FINDING EIGENVALUES 


There are versions of this program under development that will locate the eigenvalues 
and do the entire computation without user intervention. Currently, however, these versions 
are reliable only for the simpler types of profiles — usually those with only one duct — and 
are not ready to be reported. Locating eigenvalues with the standard version of the program 
is discussed here. 


The standard version of the program requires the user to find the eigenvalues. In this 
version, each time an eigenvalue is determined by iteration, the resulting value is stored and 
counted as an eigenvalue. Therefore, the user must ensure that all iterations result in good 
roots, that all required modes have been determined, and that no modes are present more 
than once. In most cases the user must expect to make more than one computer run to 
obtain this result. 


CONTROL CARDS 


The user controls the eigenvalue determination by using any of four different types 
of control cards. The first type specifies an initial value for v and an initial step size, Av. 
These are both complex numbers with a real and an imaginary part. G is then evaluated at v 
and at v + Av to start the iteration. These are essentially the Vj and Vj+1 of eq (15). If these 
two trial eigenvalues are in the vicinity of a root, the iteration will converge to that root. 


The second type of card specifies a line segment in the complex plane, along which a 
search for eigenvalues, v, is made. The end points of the line are given along with the number 
of equally spaced points at which the line is to be divided. G is then evaluated at each suc- 
cessive division point along the line until a relative minimum in IG2| is found, indicating that 
a root is nearby. The iterative process is applied to find the root. The initial step size, Av, is 
first computed to bring the second evaluation at v + Av as close as possible to the true root. 
This is done by using the point which resulted in minimum iG2| and the points on either side 
of it to determine the minimum of the parabola passing through them. If v-h,v, andv+h 
are the three points at which G was evaluated, it follows that the distance from v to the mini- 
mum of the parabola 


Av = h[G(v + h) - G(w - h)] /2[2G(v) - G(v + h) - G(v -h)]. (26) 


When the iteration is complete, the eigenvalue is recorded and the program continues 
to step along in the direction of the given line, checking again for a minimum. However, the 
stepping is resumed from the newly located root rather than from the approximate location 
where the minimum was detected. With this correction in position, the designated line does 
not have to hug the curve on which the eigenvalues are located because it is corrected at each 
eigenvalue. 


This method of finding eigenvalues has proven very successful. Its main utility arises, 
though, because the eigenvalues of the trapped modes have negligible imaginary parts and the 


search can be made along the real line. In simple profiles this can often give a successful set 
of modes on the first try. Usually, only the three initial eigenvalues need to be located by 
this means because further eigenvalues can be located by extrapolation on the previous three. 
This is the function of the third type of control card. 


The third type of card specifies the number of additional modes to be determined by 
extrapolation. The starting value of each eigenvalue is determined by extrapolating from the 
three most recently determined eigenvalues to find v. The step size, Av, is chosen as 0.0001 
times the distance between the last two eigenvalues. The exact eigenvalue is then determined 
by iteration. The extrapolation is the simple parabolic form for equal steps: 


v= 3Vn - Vell =P Wi (27) 


This method of locating modes works well when the modes lie along a smooth curve, 
as usually occurs for single ducts. But this relationship does not always occur for profiles 
with multiple ducts. 


The final control card is punched by the program when requested and contains the 
correct eigenvalue to full precision. Upon encountering this card, the program does not 
iterate, but instead evaluates G for this eigenvalue and stores this value of G as the next 
eigenvalue. A deck of such completed eigenvalues can be stored, saving the expense of 
recomputing the eigenvalues for a given profile and frequency. 


ITERATION TERMINATION 


A full description of the iteration of eq (15) should include the method of termina- 
tion. The usual criterion for stopping is that G fails to become smaller. As G approaches 
minimum size, however, round-off error can act as noise so that G is no longer a predictable 
function of v. The denominator of eq (15) can then be very small by chance, resulting in a 
large value for 6;. If this happens, the next value of v, which was as near to the root as pos- 
sible, will be far away. A much better convergence criterion is that 6; has reached a minimum 
in absolute value. In the program, iteration is stopped when |52| exceeds the previous value 
by a factor of 2. However, this criterion is not applied until three iterative steps have been 
completed, to permit the process to become well established. An upper limit of 15 iterative 
steps is permitted. We have not found an improvement on the root after 15 steps. 


SOUND SPEED PROFILE 


The normal mode program requires as inputs the depth of each layer and the sound 
speed and sound speed gradient at the top of each layer. These variables are mapped into the 
dimensionless internal variables of the program by eq (7). The purpose of the sound speed 
profile processing portion of the program is to accept the profile parameters in a form con- 
venient for the user and to translate them into the required sound speeds and gradients. 


The first function of the processing program is to make the sound speed continuous 
at interfaces. This is done simply by using the sound speed at the bottom of one layer as the 
sound speed at the top of the next. It may be necessary to compute the sound speed at the 
bottom of the layer. The necessary parameters will have been given. Occasionally a 
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discontinuity in sound speed is required, as when modeling an interface between water and 
sediment. The user indicates this by specifying the sound speed at the top of the layer. If 
left blank, the program provides the sound speed necessary for continuity. 


A second function of the processing program is to permit a layer to be defined by the 
sound speed at top and bottom of the layer rather than by one sound speed and one gradient. 
Note that the profile form as given by eq (5) is a two-parameter curve. 


The last layer extends to infinite depth, so a gradient must be specified at the top of 
it. However, this gradient can be specified by giving a depth and sound speed point below the 
last layer. The program handles this by checking to see if the gradient of the last given layer 
is unspecified. If it is, the number of layers is reduced by one, which causes the last layer to 
be only the required extra point determining the final gradient. This final gradient must 
always be negative, as is required by the boundary conditions. The program user must ensure 
that this gradient is negative and that no gradient is zero. A zero gradient will appear in the 
denominator of eq (7). 


These functions of the profile processing program are relatively simple, but an addi- 
tional capability used to model sediment bottoms greatly increases the complexity of the 
program. The capability required is to specify the absorption in a layer by adding an imagi- 
nary part to the sound speed. In older versions of this normal mode program an imaginary 
part, expressed as an absorption coefficient, could be added to the sound speed at the top of 
the layer. This imaginary part is small compared to the real part. Since tne gradient was 
assumed real at the top of the layer, the imaginary part was initially not changing with depth 
and it usually changed only a minor amount through the depth of the layer. However, this 
small change could not always be relied upon. Also Hamilton (ref 8) has published data on 
absorption gradients in sediment layers, so more precise control of this part of the sound 
speed function is needed to model sediment layers. Therefore, a more comprehensive profile 
processing routine has been incorporated in the normal mode program. This curve-fitting 
process is described below. 


The following quantities can be input for each layer depth starting at the surface: 


Depth of top of the layer 
Sound speed at top of layer 
Sound speed at bottom of layer 
Real part of sound speed gradient at top of layer 
Attenuation in loss per km at the top of the layer 
A similar attenuation at the bottom of the layer 
Density in the layer 
The density is a constant in the layer and as such requires no further curve fitting. Redun- 


dant parameters are left blank on input cards. In some cases negative values serve as flags to 
indicate specific treatment. For instance a negative value of absorption at the top of a layer 


8. Sound Attenuation as a Function of Depth in the Sea Floor, by EL Hamilton; J Acoust Soc Am, vol 59, 
p 528-535, March 1976. 


directs the program to use the same imaginary part of sound speed as occurred at the bottom 
of the previous layer. Similar flags at the bottom of a layer are discussed later. 


Absorption per Hz is given in units of decibels per km (or kiloyard). The quotient of 
absorption over frequency is used because Hamilton (ref 8) usually considers absorption (or 
attenuation) as proportional to frequency with a coefficient k. We use the symbol h instead. 
That is, 


a = hf. 


We interpret a to be in units of dB per km and f in Hz, whereas Hamilton uses dB per m and 
kHz; but the coefficients h and k remain equal. 


The complex wave number in layer i is represented as 


= wReC|C|-2 - iwlmCIC|-2. (28) 
A plane wave will be attenuated a dB per km if 


Imk; = - a/(20 000 log e) 


=-mAf, (29) 
where 
A =h/(20 000 z log e). 
By equating the imaginary part of kj in eq (28) and (29), the imaginary part of C; is found to 


be as follows: 


Y; 
ImC; = 1/A - [1/A? - (ReC)2] ~. (30) 
If a is zero, which is the case usually used in water layers, eq (30) cannot be used; but the 
imaginary part of C is then simply zero. These two cases are treated separately in the program. 


When sound speed is given at the top and bottom of layer i, the imaginary parts of 
the sound speeds are determined by eq (30) and the only curve fitting task is to determine the 
gradient y;. Solving eq (5) for ¥;, 


2 2 2, 
1 = CC, -¢ \/2Ci. (Zi4] = fis) (31) 
The gradient is a complex number since the C’s here are complex. The z’s are real. 


A second version of this computation arises if the gradient is required to be a real 
number. In this case, which is used to match older versions of the program, an additional 
parameter must be left unspecified and this parameter is Im Cj+1. This is equivalent to having 
the sound absorption at the bottom of the layer unspecified. Therefore, a negative number 
input for this parameter is used as a flag to call for this particular fitting procedure. 
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For this situation, given Re Ci, Im Cc; Re Ca and making Ne real, the determination 
of y; and Im Cj+ is not simple. When 4; is eliminated from the real and imaginary parts of 
eq (31), a quartic equation in Im Cj+] results. Rather than derive an algebraic solution to 
this equation, it is solved by iteration under Newton’s method. A good first guess at the 


solution is Im Cj+; =Im C;. Four iterations usually give an accurate root. The equation is 
Im C; (Im C,44)4 + Exes + 2(Re C41)? Im ci] (im €,,,,)7 
+ 2Re C41 Re(C))? Im Cy 
+ Im C(ReC;41)4 - (Re Ci44)2 Im(C>) = fm C; 4) (32) 
i i+] i+] i sgltelids 


The root is then found: 


The gradient, y, is next given by the relationship 
np Sy linn (Ce [ire Ce - (Im C41)? | + 2ReC; ReCi,, ImCj4, - inch] 


Because the root of eq (32) may not be exact, Im y; may not be exactly zero. This slight 
error can be transferred to Cj+ 1 by using the computed real y; to recompute C;+1. This is 
done in the program by transferring to a portion of the program already designed to do this. 


When sound speed and gradient at the top of the layer are given, the parameters 
required by the program are all given. The sound speed at the bottom of the layer is rou- 
tinely computed, however, because it may be required to make the next layer continuous. 
Equation (5) is used to determine the sound speed at depth z;+], which is the depth of the 
bottom of the layer. This is straightforward, but several complications arise. Only the real 
part of the gradient at the top of the layer is used as an input because situations have not 
arisen that require that the imaginary part of the gradient be specified. Often the attenua- 
tion is given at both top and bottom of the layer. That is, Re C;, Im Cj and Re ¥j; are given, 
plus a relationship between Re C,+1 and Im Cj+ ,. The imaginary part of the gradient, Im ¥j;, 
must be determined as well as both real and imaginary parts of the sound speed at the layer 
bottom. The derivation of this case is not trivial. 


One relationship between the real and imaginary parts of the sound speed is given by 
eq (28) and (29). From these equations at Gi we derive 


A(T = i) = 2/Ci41, (34) 
where 


T=Re C44 /Im Ci] b 


Substituting this expression for C,4] into eq (31) and equating real parts gives a quadratic 
expression for T which has a usable root of 


1 
Re(C2)T = -Im(C°) - [Im(C>)] es Re(C> )B i (35) 
where 
B= Re(C>) ~ 8 Re 7;(z;41 - Z)/A2 + 4 Re C,/A2. 
From eq (34), 


Re Cia = 2T/A(T? + 1) 
and (36) 


The gradient can now be evaluated by eq (31) to find its imaginary part. 


Equations (34) and (35) cannot be used if the attenuation at the bottom of the layer 
is given as zero. Therefore an alternate form must be used. This form is much simpler than 
the previous case, since C+] is real. 

3 yy 
Cin] = Re(C; ) [Re CG -2Re Vi (Zi4] = Z;)] (37) 


kin GR [Lien Cam Im(C>)/C2, | [Wary = z))—! (38) 


Finally, if the special case, 7; real, is specified by inputting a negative value for 
absorption, eq (31) can be used directly to give 


2 


Ci = CIC; = 27j(Zi+] = a) 6 (39) 


To evaluate the square root, let 


ie A ies ‘ 
Ci =at+t bi. 
Then 
y, Y, 
Re Cia) = (tema 1/3 (40) 
and 


NUMERICAL BREAKDOWN 


A situation arises frequently in which a very small depth function must be computed 
from the difference of two large numbers. A wrong answer results if this accuracy loss 
exceeds the word size of the computer. The best way that has been found to avoid this is to 
check for it within the program and arbitrarily replace the wrong number. In checking for 
this, a constant, called T-lim in the computer program, is compared to the argument of the 


modified Hankel functions or to the argument of the exponential function within modified 
Hankel functions. A T-lim value of 25.0 is used in the program, but a smaller number occa- 
sionally is required. The program user can alter T-lim by appropriate input cards (Key 8 = 1 
followed by a new value of T-lim). The next few paragraphs demonstrate the symptoms of 
this problem, so as to assist a user in recognizing the problem. The remainder of this section 
describes the modifications that have been made to the computer program to correct this 
loss of accuracy. 


The solid line of figure 1 shows a simple surface duct and the phase velocities of the 
first three modes at 3 kHz. For this profile, the depth function of mode | is shown in fig- 
ure 2. The solid line is the depth function as computed by a program that does not correct 
for numerical breakdown. The dashed line shows the correct depth function below a depth 
of 71 m. This result was determined from Airy functions, not from the program. Between 
depths of 71 to 100 m, the program cannot compute the depth function accurately. In the 
second layer, which starts at a depth of 100 m, the function can be computed accurately but 
it is incorrectly placed by the boundary condition that requires the depth functions to be 
continuous at interfaces. The slope of the depth function was correctly computed as indi- 
cated by the identical shape of the three depth functions in the second layer. The shape is 
such as to make the correct depth function continuous in slope across the interface. 


The breakdown in accuracy at a depth of 71 m occurred when ¢ had a value of -8.4. 
(¢ is given by eq (7) and is the argument of the modified Hankel functions.) A negative value 
of ¢ occurs when the mode phase velocity is less than the speed of sound. Since the ray of 
the same phase velocity cannot reach such a region, the sound field there is a diffracted field. 
The mode depth function is therefore small at such depths. In the figure, the depth func- 
tion amplitude at the breakdown point is about 7 orders of magnitude (or in terms of propa- 
gation loss, 140 dB) down from its maximum. Equations (62), (66), and (68), which will be 
given for the modified Hankel functions, indicate that the argument of the exponential term 
is 2/3(8.4)3/2, or 16.2. The functions hy and h9 will thus be about 10/7 in magnitude at a 
depth of 71 m. These large values and their small difference account for the approximate 
accuracy loss of 14 decimal places, which is the general accuracy of the modified Hankel 
functions. 


Incorrect behavior in the depth function usually occurs when ¢ is about -8.4. In 
some more complicated profiles, however, where accuracy is lost in row reduction of the 
determinant, the depth functions may become incorrect at values of ¢ that are less in abso- 
lute value. When this problem occurs it can be diagnosed by plotting the depth function of 
the mode and noting the steep positive slope through some depth interval as in figure 2. 
When that occurs, the value of T-lim should be decreased. 


Incorrect depth functions can cause errors in propagation loss computations in two 
ways. In figure 2, the solid-line depth function, because of its large size, can cause losses to 
be too low at a depth of around 100 m. The second error would occur if the duct were 
deeper, say 110 m. At this depth the erroneous segment of depth function in figure 2 would 
reach a value of about 107!, where it would be larger than the correct lobe of the depth 
function near the surface. With this extra area under the curve, the normalizing factor would 
be increased significantly and would reduce the size of this entire depth function. Thus, 
losses near the surface would be larger because of the loss in size of mode 1. 
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Figure 1. A two-layer sound speed profile for a 
surface duct. The phase velocities of the first 
three modes at 3 kHz are marked. The broken 
line shows a modification of the upper layer to 
prevent numerical breakdown in mode 1. 
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Figure 2. Depth function of mode 1 at 3 kHz, showing error in the computed 
function. The true function cannot be computed without increasing computer 
word length, but the corrected value can and it will not cause a large error in 
the mode sum. 
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The standard correction to mode | is shown in figure 2 by the dot-dashed line. In the 
depth interval where § <-8.4, the function is set to zero. The values of the depth function 
at greater depths result from a modification in the values used in the determinant. 


The corrected values in figure 2 are not equal to the true value of the depth function, 
but they are small enough that they do not alter the propagation loss to a tenth of a decibel 
when a full set of modes is used. When the source or receiver is at a depth where such cor- 
rections are necessary, the mode can be omitted from the computation. Thus, properly 
omitting modes would solve the above problems except for the cases where the normalizing 
factor, D, is affected. In these cases, losses cannot be computed accurately without the 
corrections. 


PROGRAM MODIFICATIONS 


The modification is approximately equivalent to modifying the sound speed profile 
as shown by the broken line in figure 1. In effect, the sound speed is not allowed to become 
enough greater than the phase velocity of the mode being considered to cause problems. 


The limitation on ¢ is accomplished at three different places in the normal mode pro- 
gram. It is not clear that this is the best way to handle the problem and it may be redundant, 
but it appears to be an adequate solution. These three corrections will be described next. 
Finally a correction to the determinant program is described which is necessary because the 
limiting of € can cause false zeroes in the determinant. 


In the subroutine SETUP the elements of the determinant are computed by deter- 
mining ¢ at the top and bottom of each layer and then calling the modified Hankel function 
program. At the top of each layer, Re §¢ is set to -7.5 if its value was less. However, this is 
done in an iterative loop in which the real part of w/C; in eq (7) takes on a sufficiently larger 
value while its imaginary part is fixed. This is done to retain the absorptive properties of a 
layer when its sound speed is in effect being reduced. It has been found unnecessary to make 
the above constant, -7.5, a function of T-lim which the user can vary, because an oversized 
value at the top of a layer is not as critical as at the bottom. At the bottom of a layer, several 
tests are made. If the real part of ¢ has decreased past the limit at some depth between the 
top and bottom of the layer, it is set at the limit. This limit, called S-lim in the program, is 
related to T-lim by the relationship 


s=-(1)2/3 (42) 


where S and T are the two limits. If Re ¢ is less than -7.5 throughout the layer, it is simply 
set at -7.5. Such a layer has negligible effect on a mode. 


In program MAIN at the location where depth functions are computed for given 
depths, a process similar to that above is used. To evaluate the depth function in a given 
layer, ¢ is first evaluated at the top of the layer. The real part of § is then limited as in the 
program SETUP above. Next § is evaluated at the given depth by adding the depth-dependent 
part onto the value at the top of the layer which may be the modified value. If this final 
value is less than S-lim, the depth function is set to zero. If it is greater than S-lim, the func- 
tion is computed in the usual way. 
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The imaginary part of ¢ can be large if the eigenvalue has a large imaginary part or if 
the speed of sound in the layer has a large imaginary part. When this happens the imaginary 
part of 2/3 ¢3/2, which appears as an exponential in the modified Hankel functions, may 
become large in absolute value even though Re § has been limited. A final check is therefore 
made before the exponential is computed. If Im 3/2 is greater than T-lim, ¢ is reduced in 
amplitude to the size at which it will equal T-lim. The angle of ¢ in the complex plane is 
preserved. 


This limitation of the exponential can be viewed in another way. In a following sec- 
tion the two components of the modified Hankel functions, F; and F, eq (68) and (69), 
have exponential terms whose arguments are equal and opposite in sign. When these argu- 
ments have magnitude of 2/3 T-lim, they differ in size by 15 decimal places, which is near 
the 18-decimal-place word size of the machine. The ability to compute the difference in 
these two terms is essentially the same as the ability to compute the depth function accurately. 


PREVENTING ZEROES IN THE DETERMINANT 


Placing limits on § can cause problems in the determinant because § may be set equal 
to S-lim at several interfaces. The equations that arise for matching boundary conditions 
may then be identical for these interfaces and may therefore fail to be linearly independent. 
The triangularized determinant will thus have zeroes on the diagonal at positions equivalent 
to interfaces that do not have real physical importance for the mode. These will prevent 
location of the significant “‘zeroes” or roots. These artificial zeroes must be removed. 


The artificial zeroes are detected and removed in the subroutine DET, which evaluates 
the determinant. If four elements from the matrix have the configuration 


a b 
Cc d 


and c is to be set to zero by row reduction, d will be replaced by a value, x, as follows: 
x =d -be/a. 


If d is located on the diagonal, complete loss of accuracy is checked for by computing 


s= Ix21/ 1421. 


If s is less than 10734, x is not used; instead, d is replaced by 107!74. Note that this substi- 
tution will occur when x is zero, thus preventing zeroes on the diagonal. The power of ten, 
-17, is chosen to be near the total word size of 18 decimal places. 


The above substitution prevents sudden jumps in the value of the determinant when 
all precision is lost at one step in the evaluation. This is important for the mode search rou- 
tine which detects roots by looking for minima in a series of values of the determinant while 
one parameter is incremented slowly. A sudden jump will often produce a relative minimum 
which will be falsely interpreted asa root. At true roots, one or more elements along the 
diagonal are small, but not as small as those checked for here. 
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REFLECTION COEFFICIENTS AND OTHER AUXILIARY OUTPUTS 


Once the depth functions of a mode have been determined, it is relatively easy to 
compute reflection coefficients at any interface. Therefore, a subroutine called RCOEF has 
been added to the program which will compute and print out reflection coefficients if 
requested by the use of control key 3. If key 3 is set to 1, the reflection coefficients at all 
interfaces are computed. If set to a number, n, greater than 1, the coefficient is computed 
at the nth interface only, where the surface is the first interface. 


The printout includes the phase as well as the amplitude of the reflection coefficient 
and the grazing angle. The grazing angle, 0, of the equivalent rays is computed from the 
mode phase velocity and the sound speed, c, at the bottom of the layer, by Snell’s law: 


6 = cos”! (c/v). 


The grazing angle is computed only if the phase velocity is greater than the sound speed at 
the interface, since otherwise the equivalent ray does not reach the interface. 


The reflection coefficient is derived, following Bucker (ref 9), by assuming that an 
isospeed layer exists for a small depth just above the interface. In this layer the depth func- 
tion can be written as 


f(z) = AcllZ + Bellz, (43) 
where 1, the vertical component of the mode wave number, is given for mode n by 
2D) 
eke Lye (44) 
in a n 
and 
SCT 


where Cpj is the sound speed at the bottom of layer i. The derivation now consists of identi- 
fying A and B as the pressures of the upgoing and downgoing waves at the bottom of the 
layer; thus the reflection coefficient 


R=A/B. 


A and B are evaluated by making f and its derivative at the interface between this small iso- 
speed layer and the regular profile continuous with the normal mode depth functions. The 
thickness of the isospeed layer is then allowed to approach zero, giving the desired value of 
R. If Fand F’ are the normal mode function and its depth derivative at the interface depth 
defined by eq (9) and (19), the reflection coefficient resulting from the above derivation is 
as follows: 


R = GilF + F’)/(IF - F’). (45) 


This coefficient is a complex number. Loss per reflection is given by 20 times the 
log of the absolute value. The phase gives the phase shift that an equivalent ray would 


9. Sound Propagation in a Channel with Lossy Boundaries, by HP Bucker; J Acoust Soc Am, vol 48, 
p 1187-1194, November 1970. 
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experience upon reflection. Figure 3 is an example of the use of this computation. It shows 
phase and amplitude of the reflection coefficient in shallow water over a sandy-silt sediment 
lying over rock. The frequency is 1500 Hz. Reflections are given only at discrete points 
determined by the individual modes. 


The model in figure 3 is for a liquid bottom. That is, no rigidity is supplied in this 
program and the sound speed, density, and attenuation determine the reflection coefficients. 


The reflection coefficients computed by eq (45) can be closely approximated by 
dividing the mode attenuation by the loop length of the corresponding ray. The loop length 
must be determined from ray theory for the ray of the same phase velocity or vertexing veloc- 
ity. However, an interesting analog of the ray loop length is the intermode interference length. 
This is discussed by Guthrie (ref 10). Specifically, if the difference between eigenvalues, Re Nip 
for two adjacent modes is AX, the interference length 1 = 27/AX. This distance will usually 
equal the ray loop length for some ray with phase velocity between that of the two modes. 


As each mode after the first is computed, the length, 1, is computed and printed out. 
Also routinely printed out for each mode is the mode damping or mode attenuation coeffi- 
cient, in units of dB per km. This attenuation, aj, is computed from the relationship 


a; = -1000 Im A; log; ge 
= — 8686 Im );. 
This quantity multiplied by range gives the damping of mode i, in dB. 


10. The Connection Between Normal Modes and Rays in Underwater Acoustics, by KM Guthrie; J of Sound 
and Vibration, vol 32, no 2, p 289-293, 1974. 


SOUND SPEED, m/s 


read 1530 1540 1550 2.5 250 
fea} 

ke = ui 
2 \o fe 
5 2.0 200 © 
100 5 w 
a) 
uu . 
S 1ST SS RUMEN 1.5 150 F 
. Ww 4, — 
= 200 0.3 m THICK ic ASe 55 
ui 2ND SEDIMENT ti 1.0 100 
2) LAYER e < 
300 m THICK D a 

300 S05 50 

(0) 0 


FINAL LAYER 
400 0 2 4 6 8 10 


GRAZING ANGLE, DEGREES 


Figure 3. A shallow-water profile with resulting phase and amplitude of the reflection coefficient at 1.5 kHz. 
Parameters at the top of the sediment layers are as follows: 1st layer — c = 1606.45 m/s, y = 1.5371, a= 
0.18 dB/m, p = 1.68; 2nd layer — c = 1684.0 m/s, y = 1.5s7!, a = 1.10 dB/m, p = 1.91; final layer — y = -0.1. 
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COMPUTATION OF THE MODIFIED HANKEL FUNCTIONS 


Most of the computer time required to determine eigenvalues and compute depth 
functions is spent in evaluating the modified Hankel functions of order 1/3. For this reason, 
minimizing computer time in evaluating these functions is desirable. Gaining as many places 
of accuracy as possible is even more important. The average normal mode computation will 
have many modes that can be determined to far greater accuracy than is required to obtain 
0.1 dB accuracy in the propagation loss. However, there are usually some and often many 
modes in which many places of accuracy are lost in evaluating the determinant. Therefore, 
maximum accuracy in the modified Hankel functions is required to extend the range of 
cases for which computations can be carried out successfully. 


Optimization of the program is a function of the computer word length. The pro- 
gram given in this report is for the UNIVAC 1110 with 60 bits word length in double preci- 
sion or 18.1 decimal places. This section gives the equations and computational techniques 
that are required to optimize this program for different computer word lengths. Complete 
details of the functions are given in reference 11. 


The Airy functions Ai(Z) and Bi(Z) can be used instead of the modified Hankel 
functions hy and hj. However, since ho is ideally suited to matching the boundary condi- 
tions at great depth as formulated in this normal mode program, hy and hp are used here. 
The relationship between them is as follows: 


hy @) = k [Ar (@z)=1 Bi €z)) (46) 


k* [Ai (-z) +i Bi (-z)] , (47) 


hy (Zz) 
where 


k = (3/2)2/3 (1 - i,/3/3), and k™ is the complex conjugate of k. 


In this section z will be the argument of the functions hy and hy. For small values 
of |z|, h and hz are computed by power series expansions. For large values, an asymptotic 
expansion is used. In the past the asymptotic series was expanded directly. However, a 
continued fraction expansion has been found to give both shorter running time and better 
accuracy. 


Figure 4 shows a line in the complex plane which divides the plane into two parts. 
For values of z within the line, the power series method is used. When z is outside the line, 
the continued fraction method is used. This line is a function of computer word length, and 
the method of determining it will be given after the two methods have been treated. The 
accuracy of the methods is also treated. 


The program has a parameter called IH in the FORTRAN call statement which con- 
trols which functions are computed. If IH is set to zero, both functions and their derivatives 
are computed. If IH is set to 1, only the functions are computed. If it is set to 2, only hj 
and its derivative are computed. 


11. Tables of the Modified Hankel Functions of Order One-Third and their Derivatives, Harvard University 
Computation Laboratory; Harvard University Press, Cambridge MA, 1945. 
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Figure 4. Line in complex plane dividing the arguments for which the 
modified Hankel functions are computed by (1) power series (inside) 
and (2) asymptotic expansion evaluated by continued fractions 
(outside). 


POWER SERIES EXPANSION 


In this expansion hy and ho are given by two auxiliary functions f and g as 


hy @) = g+i@G)!/2@-2f (48) 


hy (2) = g-i@3y !/2 (g-2f) . (49) 


The auxiliary functions are given by the expressions 


M 
ve a x (50) 
m=0 
M 
3 = Ee Bey ee (51) 
m=0 


where X = 23, A = 21/3/[17(2/3)] and B = 2!/3/[32/3 Pr (4/3)]. The derivatives h} (z) and 
h» (z) can be derived by straightforward differentiation of eq (50) and (51) to give 
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ff = -Az2 » exam (52) 
m=0 
M 
Pre) i chee t (53) 
m=0 


The coefficients of eq (50) - (53) are given by recursion relations where ag = Il, a = 1/3!, 
Ay Giarll 9 cil me Sa ohie.di.o 1/0) 


am = —a-1/(3m) 3m-1) (54) 
bg = 1, by =-2/4!, by =+2 - 5/7!,b3 =-2 - 5 - 8/10! 


Dra Oren omm) Gmail) (55) 


m 
cg = 3/3!,c) =-6- 1° 4/6! 


GA HC. 7d Gan) (56) 


m 
doz lid] =-4- 2/4! 


diy = -dy)-1/3m (m-2) (57) 


m 


It is important for efficient computation that the number of terms M be no larger 
than necessary. In the current program the same value of M is used in all four sums. This is 
done because the optimum number never differs by more than one in the four cases and the 
determination by table look-up of four M’s often would take longer than computing any 
unnecessary terms. M for each series is determined so that adding additional terms will not 
change the answer. Then the most stringent of the four conditions is tabulated and used. 


A precise determination of the number of terms to use requires a knowledge of the 
size of the largest single term in the sum. When a term is smaller than this by a factor which 
is the power of 10 equal to the number of decimal places in the computer word size, it 
cannot affect the sum. We ignore the fact that a sum of small terms might be significant. 
This, then, defines the truncation point. Let m be the number of the largest term in the 
sum, k the number of terms to be used, and h the number of decimal digits in the machine 
word. Then for a given k, the largest absolute value of the argument z that can be used to 
compute g’ is given as 


m k 
eo. = el Gee 10 . (58) 


The power of ten can be replaced by 2 raised to a power of the number of binary bits in the 
computer word if preferred. The coefficient d of eq (57) is used. Each of the other three 
should also be tried, to find the smallest number of the group for a given k. Equation (58) 
can be solved for |z|, giving 


log |z| = (log dy - log d,, + h)/(3m - 3k) . (59) 
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A simple computer program given in appendix C will find |z| for each value of k from | up 
to the maximum number of terms desired. The largest term, m, is easily determined because 
from one k to the next m will remain the same or increase by 1, so it is only necessary at 
each step to check term m+1 to see if it is larger than term m. 


The FORTRAN subroutine HANKEL given in appendix A uses the above power 
series method to compute hy and hj for small arguments. The coefficients a, b, c, and d are 
given in lists by that name. The truncation points are given in the list called ZMLA2, which 
lists values of |z|2 determined by eq (59) or the three similar equations. 


ASYMPTOTIC SERIES EXPANSION USING CONTINUED FRACTIONS 


When the argument z falls outside the curve in figure 4, hy and hj can be computed 
more efficiently or more accurately by asymptotic series than by power series methods. 
Reference 11 gives information on branch cuts and regions of validity of the two forms of 
the asymptotic solution (Stokes’ phenomenon). Here we will give computing formulas that 
comply with these requirements, without discussing them further. 


Since a given expansion is valid in one or more quadrants, we choose complete quad- 
rants as regions. For z in quadrants 1, 3, or 4 use 


hy (z) ~ exp (Sm i/12) Fo (z) (60) 
h5 (z) ~ exp (-7 i/12) G5 (z) (61) 


For z in quadrant 2 use 


hy (z) ~ exp (57 i/12) F> (z) + exp (117 i/12) 7, @) (62) 
h5 (z) ~ exp (-7 i/12) G5 (z) + exp (-777 i/12) G, (z) (63) 


For z in quadrants 1, 2, or 4 use 


hy (2) ~ exp (-5m i/12) Fy (2) (64) 

hi, (2) ~ exp (m i/12) Gy (2) (65) 
For z in quadrant 3 use 

hy (z) ~ exp (-5m i/12) Fy (z) + exp (-117i/12) Fy (2) (66) 

hi, (z) ~ exp (m i/12) Gy (z) + exp (77 i/12) G5 @) (67) 


The four auxiliary functions follow: 


M 
Fy (2) = K27!/4 exp (2: 23/2/3) » Cy (68) 
m=0 
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M 
Fy @) = kz "4 exp (-2i73/2/3) Sc, ¥™ (69) 
m=0 
M 
eZee (ilz/2/3) » Din X™ (70) 
m=0 
M 
CHG) zy exci) De ae (71) 
m=0 


Gy @) 


3/2 


where X and Y equal + iz respectively, and 


K = 2!/3 31/6 1/2 = 9.853 667 218838951 


The coefficients C,, and D,, are again computed by recursion relations where Cg = Do = 1: 
Cm = Cm [9 (2m-1)? - 41/48m (72) 

and 

B=) Dl 9(2m=1)~ 116) 48min (73) 


Square roots of z are to be taken so that the real part of the root is always positive and the 
imaginary part has the same sign as the imaginary part of z. This applies also to fourth 
roots. The three-halves power is obtained as the product of z and its square root. 


The summations in eq (68) - (71) can be done as indicated or evaluated by contin- 
ued fractions. When done as indicated they are asymptotic series, and care must be taken to 
truncate them at the term of smallest magnitude, if this term is reached, because adding 
more terms will reduce the accuracy. Since the largest term in these series will always be 1, 
the series can be truncated if the terms become less than 107) in magnitude, where h is the 
number of decimal digits in the computer word. 


Continued Fraction Expansion 


The method of continued fractions is more effective in evaluating these asymptotic 
series, and it is used in subroutine Hankel in the FORTRAN program in this report. The 
coefficients are stored in lists entitled C4, C5, D4, and DS. In the remainder of this section 
the continued fraction technique is presented, along with the method of determining 
coefficients. 


The continued fraction has the form 
F (x) = bo + ay 
Xoo by + a9 : (74) 
XGct 1) Ts aoe 


It is to be used to evaluate a polynomial 
M 
Pe) = » CeO (75) 
=0 


This polynomial can represent any of eq (68) -(71). One of three standard forms for con- 
tinued fractions, this form is used because it has two coefficients at each stage and therefore 
is equivalent to an asymptotic series of twice as many terms. This reduces by half the num- 
ber of divisions required. Since complex divisions are lengthy, requiring six real multiplica- 
tions and two divisions, this is the only standard form of the continued fraction that can 
compete in computer time with the asymptotic series. 


The coefficients a; and bj of eq (74) must be determined from the coefficients Cy). 
The usual technique is to express P as a rational function, then use the continued fraction to 
evaluate the rational function. The determination of the coefficients can be done in these 
two steps or by a second method which goes directly from power series to continued frac- 
tion coefficients. The second method is preferable because the loss of accuracy is more in 
the first. But since the first method is more easily understood, each method will be given; a 
computer program is included in appendix C which will determine coefficients by the sec- 
ond method. 


Let M in eq (75) be an even number so that 2N=M. (An additional unnecessary 
term of the series can always be used.) The rational function will have the form 


N N 
R(x) © k » e x! » oe (76) 
i=0 i=0 


where eg = fg = 1 and k = Cg. The coefficients ej and fj are evaluated from a set of linear 
equations which can be described by displaying a particular case. For N = 3 they are as 
follows: 


aN Oven OO 

OS 0 (Gc Oo C5 

© 0 1 G&G G _ | ve 
0 0 0 Ee G Ch 

OVO) “Cre, oC Cs 

OQ ONO Cy KeANC. C6 


With e; and fj thus determined, R(x) is equivalent to P(x) through the first M + 1 terms. 
R(x) can now be evaluated exactly (except for round-off error) using a continued fraction 
of the form F (x) of eq (74). 


Di 


Rather than R(x), however, a similar expression in y = 1/x is the form that is well 
suited to evaluating asymptotic series. This expression is obtained by dividing each term of 
R(x) by XN_ The order of the coefficients can now be reversed and a simple algebraic oper- 
ation can yield a value of | for each of the two initial coefficients and a new value for k. We 
will call this new rational function with renamed coefficients R(y). It will have the form of 
eq (76) but different coefficients, say e and f instead of € and f. 


The coefficients a; and b; are determined from e; and fj; by a recursive formula which 
involves constructing an n X n triangular matrix Q with elements qi,j as follows: 


Lo) FEO 

dy i = (-e9 fp/ay eon linea opera N eo 
where Cy = 1, giving a, and 

Dy ah) os 
The second row: 

Go = Cho san oy qy p/ag is 3 De Nid 
where Go 1, giving a7, and 

ly) Gin oy 3) 


Elements outside the matrix are assigned a value of zero. The remaining rows for m = 3 to 
N are as follows: 


dm,i = (m-2,i7 Im-1,i+1 ~ &m-1 Im-1,i)/4m i= m,mtl,...,N, 


where diy m = |, giving ap, and 


bm = Im-1,m74m,m+1 


The second method determines the continued fraction coefficients a; and bj; directly 
from the asymptotic series coefficients cj. This method is preferable to the first because the 
loss of accuracy in inverting the matrix in eq (77) can be more than the loss in this second 
method. 


It has been pointed out™ that the second method is probably a variant of the 
Viskovakoff algorithm described by Khovanskii (ref 12) and as such is unstable — subject to 
accumulation of errors. However, it is sufficiently stable to obtain the required coefficients. 


*Private communication with AN Stokes, CSIRO, Wembley, Western Australia. 
12. The Application of Continued Fractions and their Generalizations to Problems of Approximate Analy- 
sis, by AN Khovanskii; a monograph in Russian, 1956. 
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The coefficients are derived as follows. The well-known recursive relations that give 
the Nth stage of a continued fraction as a rational function are used (ref 13). 


Fy(y) = Ay(y)/Byy) « (78) 
where 
N 
An = » ej y! 
i=0 
= (y + by) An-|] + an An-2 (79) 
N 
By S » i y! 
i=0 
= (y + by) By_} + an By-p > (80) 


in which A_; = 1, Ag = bo, B_j = 0, and Bg = 1. Again y = 1/x. The long division indicated 
in eq (78) is then carried out, giving a quotient in terms of aj, bj, and y that can be equated, 
term by term, to the first 2N-1 terms of the asymptotic series. 


The long division is carried out with Ay and By written in descending powers of y. 
The quotient is then in descending powers of y or ascending powers of x. Fortunately, the 
first 2N+1 terms determined for any N are identical to the same initial terms for any larger 
value of N. This will be proven later. The first few equations obtained from the division are 
as follows: 


ag = Gh 
ay = (Cy 

ay Og — & 

ay (>;-*2) = (Ca 

ay (2 ay by -b3 +29ba) =C) (81) 


From these equations a; and b; can be determined, since the coefficients Cj are 
Known. However, a simpler method is available. 


The long division indicated in eq (78) can be carried to 2N+1 valid places; but be- 
yond N+1 places, terms from the original dividend are no longer entering the remainder. 
Therefore terms in the later part of the quotient have a simplified form. Since term n+1 of 


13. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, ed by M 
Abramowitz and IA Stegun; National Bureau of Standards Applied Mathematics Series, vol 55, p 19, 
1964. 
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the quotient is equal to the nth asymptotic coefficient, designate it C,,. Note that the C’s 
are numbered from 0 to N. Let the coefficient of y™ in By be By m. Then 


N 
CNeene » BN ,N-i CN4j-i ° SJ SIN (82) 
i=l 


Here the C’s are numerical constants. The unknowns, the a’s and b’s, are in the terms of B. 
Suppose that these unknowns have been determined up ton = N-1. Then eq (82) will 
contain two unknowns, ay and by. By using eq (82) for} = N-1 and N, the unknowns can 
be evaluated. The index N can then be increased by 1 and the process repeated. The 
process can start with N = 2 if aj, bg, and by are provided, but these are easily determined 
from eq (81). The terms of By are determined from eq (80), which gives for each term 


B 1 +b, B +a,B (83) 


nm ~ Bu-1,m- n-l,m n-2,m ° 


Any By, m is zero if m is greater than n. 


When j = N-1 is used in eq (82) in the process described above, the coefficient of 
the (2N-1) power of x is being evaluated. This term is expected to contain ay and by, but 
— as will be proven later — because the coefficient of by is zero, ay is the only unknown in 
a linear equation and can be easily evaluated. The next term determined with j = N contains 
ay and by, but now only by is unknown and is easily evaluated. 


As an example, the Crs through n = 10 are listed in table 1. These are the asymp- 
totic series coefficients given by eq (72). The corresponding a,,’s and by’s as determined 
above are also listed. A more complete list of the a’s and b’s can be obtained from the 
FORTRAN program in appendix C. 


Table 1. Asymptotic series coefficients, C,,, and the corresponding 
continued fraction coefficients, ay and by. 


i 0 

1 0.10416 0.10416 -0.80208 
2 0.08355 -0.58764 ~2.28555 
3 0.12823 ~2.29072 -3.77864 
4 0.29185 -5.11525 -5.27462 
5 0.88163 -9.06285 -6.77193 
6 3.32141 
7 14.99576 
8 78.92301 
9 474.45154 

10 3207.49009 
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A FORTRAN program to compute the continued fraction coefficients for the series 
given by eq (72) is given in appendix C. This program can be easily modified to determine 
the other set. 


Two Proofs 


In this section proofs will be given of two facts used in the previous section. Follow- 
ing this, the number of terms required, the accuracy, and similar topics will be discussed. 
To prove that the first 2N+1 terms of the quotient Ax/ByN are equal to the same terms 
when N is a larger integer, use long division on eq (79) and (80) to obtain 


Ay/By = Ay-1/By-] + an (An-2 BN-1 - An-1 BN-2)/(@n By_p) - (84) 


If the first quotient on the right is to have terms equal to the quotient on the left up 
through term 2N-1, the remainder must have no terms with y to a higher power than 
-(2N-1). The final divisor, By By_1, contains y to the (2N-1) and lower powers. There- 
fore, the proof is complete if the numerator of the remainder is a constant. To show this, 
use eq (79) and (80) to evaluate BN_; and Ay]; it can be shown that 


Ay-2 By-1 - AN-1 BN-2 = -8Nn-1 (AN-3 BN-2 - AN-2 BN-3) 


(-1 yN an-] an-2 5.0.0 ay 


The right-hand product is obtained by repeatedly applying the middle result. The product 
of a’s is a constant, completing the proof. 


The second proof required is that in the quotient of AN/BN, Con (the coefficient of 
y72N) will contain no aj or bj to higher than term N and Cyn 4] (the coefficient of y72N-1) 
will involve no a; to higher than term N+1 and no b; to higher than term N. 


The first part is intuitively obvious. Since from the preceding proof Cay will be the 
same when derived from the ratio Ay/B, for any x as long as it is N or greater, we need 
consider only the case where x is N. But since from eq (79) and (80) Ay and By contain no 
a’s or b’s of greater than term N, Con cannot contain a’s or b’s of higher terms. 


By the same argument Cyy4] can contain no a’s or b’s to higher terms than N+1. 
There remains to be proven only that b+) cannot exist in CpN41 or that its coefficient, 
which we will call E, is zero. Applying eq (82) for N+1 and j = N gives 


N+1 


CN+1+N = 7 » BN+1,N+1-i CN+1+N-i - 
i=l 


E, the coefficient of b,+) in this expression from eq (83), takes the form 
N+1 


Eun » By N+1-i C2N+1-i - 
i=l 
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But by choosing j = N in eq (82) we see that terms | to N for Cyn equal terms 2 to N+1 in 
E, so 


E= ~ By N Con t+ C2N : 


However, since BN,N> the coefficient of y™ in By, is always 1 by eq (80), E=0. Therefore 
bn+ 1 does not exist in Cp N+]. 


Number of Terms 


The number of terms or stages to use in the continued fraction was arrived at by a 
trial and error process. For a given number of terms, a real positive argument was decreased 
until the accuracy began to drop. The magnitude just before this drop was considered to be 
the optimum point to increase the number of stages by one. Because the argument to the 
continued fractions is 73/2, we took the larger of the magnitudes of the real and imaginary 
parts of z3/2 as the test number. This number is then compared to the 3/2 power of the 
points determined along the real axis by trial and error. 


The above method appears to work well although it involves no thorough under- 
standing of the way complex numbers affect the successive convergents of a continued 
fraction. Table 2 shows the points down to which a given number of stages gives full accu- 
racy for positive real arguments and lists the 3/2 power of these numbers as used in the 
FORTRAN program list called ZMLAS. 


Division Lines 


The power series method is now to be used for small arguments and the continued 
fraction method for large arguments. The exact dividing line between them is needed. The 
division line of figure 4 was arrived at by computing the functions along rays from the ori- 
gin, using both power series and continued fractions. The number of decimal places to 
which the functions determined by the two methods agree tends to reach a maximum at 
some distance from the origin along each ray. At distances short of this maximum we can 
assume that the continued fraction method is less accurate than the power series. At dis- 
tances beyond the maximum, the power series is assumed to be less accurate. The maxi- 
mum therefore indicates the ideal place to change from one method to the other if the 
decision is to be based solely on accuracy. This method was used to determine figure 4. 


A complication arises, however. Along certain rays from the origin, hy and its deriva- 
tive reach a maximum number of places at very different distances from h9 and its deriva- 
tive. The principal problem is at +60° but persists from about 30° to 90°. At 60°, hy is 
small in magnitude and ho is large. The power series method cannot compute the small 
values accurately due to loss in accuracy in subtraction in eq (48). The accuracy of the 
continued fraction for h9 is poor at 60° because eq (69) becomes a nonalternating series and 
continued fraction approximations are not known to improve the accuracy of nonalternat- 
ing asymptotic series as they do for alternating series. 


A reasonable solution to this problem is to compute hy by continued fractions and 
hy by power series for arguments at these angles and magnitudes from 4 to 10. However, as 
will be shown later, the above solution has not been employed at this time since this area is 
not of great importance for normal mode computations. Instead, the argument was chosen 
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Table 2. Cut-off points for determining the number of stages in the continued fractions. 


Real Argument Program Test Value 


Number of Stages x3 
1 10? 
2 80 715.0 
3 35 207.0 
4 22 103.0 
5 13 47.0 
6 11 36.4 
7 9 27.0 
8 8 22.6 
9 7 18.5 
10 6.5 16.6 
11 6 14.7 
12 5.8 14.0 
13 Ss) 12.9 
14 58} 12.2 
15 5.1 11.5 
16 4.9 10.8 
17 4.5 9.5 
18 44 9.2 


that gave equivalent accuracy for the two methods. Along 60° this minimum accuracy is 9 
decimal places. 


The following relationship exists between hy and hp for positive and negative values 
of the imaginary part of the arguments: 


hy (2*) = [hy (@]” , 


where the * means complex conjugate. Thus, the above discussion at 60° can be translated 
to ~60°. Also, the functions actually need to be computed only in quadrants I and II. They 
could then be evaluated in quadrants III and IV by the above relationship. The above rela- 
tionship explains the symmetry of figure 4 about the real axis. 
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COMPARATIVE ACCURACY 


The accuracy of the three methods — power series, asymptotic series and continued 
fractions — has been determined on a CDC computer with 48 bits or 14.4 decimal places of 
accuracy in the floating point word. Since this differs from the double precision word 
length of 60 bits or 18.1 decimal places that applies to the preceding part of this report, 
these results are for comparative and illustrative purposes only. 


Accuracy is determined by computing the functions and either comparing the an- 
swers for the several different computing methods or computing the wronskian. The 
wronskian is a constant given by the relationship 


hy hy -hy hy = -1.45749544104i = -1961/3/n (85) 


The wronskian will determine the accuracy of the functions if it can be computed without 
loss of accuracy. If the two products in it are large, though, accuracy will be lost in the 
subtraction. This generally happens for arguments near the negative real axis. Here accura- 
cy must be determined by comparing answers from different methods. The accuracy of the 
functions and their derivatives will generally be about equal. 


Figure 5 illustrates the accuracy that is obtained in different parts of the complex 
plane of the argument, z, by using the power series method. On the inner contour, the 
functions h, and hp and their derivatives have 12 places of accuracy. On the outer contour, 
the accuracy is 11 places. As expected, the accuracy is best for arguments of small magni- 
tude. The accuracy remains best in directions from the origin in which the functions are 
large in magnitude. This is because less accuracy is lost in subtraction. Accuracy must be 
lost when individual terms of the series are large but the sum is small. 


Figure 6 shows accuracy contours for the asymptotic expansion with both the direct 
and continued fraction evaluation of the series. Here, the best accuracy is obtained for large 
arguments, and accuracy decreases toward the origin. As can be seen, each of the two meth- 
ods is better in some directions from the origin. The choice of methods then depends upon 
which directions are of most value to the normal mode program. The dots on the figure 
show the locations at which the functions were evaluated in a typical surface duct run. 
Although arguments can lie anywhere in the plane, most of them follow this pattern. They 
lie just above the negative real axis and in a narrow angle above the positive real axis. The 
continued fraction method is distinctly better on this positive side. Since computing time 
also favors the continued fraction method, it is clearly the method to use. 


If the 12-place accuracy contour from figure 6 lies inside that for figure 5 at some 
angle from the origin, 12 places can be obtained at any range along this angle by using either 
power series or asymptotic expansion in the interval of overlap. If the asymptotic expan- 
sion contour lies outside the other, there is an interval in which 12 places cannot be ob- 
tained. Only some lesser number of places can be obtained in this interval. These contours 
apply when both functions and their derivatives are all computed by a single method. As 
mentioned earlier, increased accuracy could be obtained in some areas by computing the 
two functions by different methods. 
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11 PLACE ACCURACY 


12 PLACE ACCURACY 


Figure 5. Locus of arguments for which the power series evaluation of the 
modified Hankel functions gives 12 and 11 decimal places of accuracy for 
a computer word length of 14.4 decimal places. 
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Figure 6. Locus of arguments for which the direct and continued fraction 
evaluation of the asymptotic series gives 11 and 12 decimal place accuracy. 
The arguments at which the modified Hankel functions were evaluated in a 
typical normal-mode run are shown. 
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MODE FOLLOWER PROGRAM 


Appendix D lists the Mode Follower Program in FORTRAN. It is not a part of the 
general normal mode program, but is related in that it uses some parts of the general pro- 
gram. The purpose of the mode follower is to trace a given eigenvalue as some parameter is 
varied. This parameter is usually frequency, but any profile parameter can also be varied. 
The eigenvalues at a given set of parameters are discrete points. By permitting the parame- 
ter to vary, the eigenvalues become a set of lines, and this often clarifies their behavior at the 
fixed points. Figures 7-9 illustrate this. 


Figure 7 is a sound speed profile consisting of two ducts. Figures 8 and 9 show the 
real and imaginary parts of some eigenvalues of the profile over a range of frequencies. The 
imaginary parts are expressed as mode attenuations. The figures show a region where both 
ducts are exerting an influence on the eigenvalues. The broken lines show the location of 
eigenvalues for a profile that consists of only the upper duct of figure 7. Considerable time 
could be spent studying the interaction between the two ducts, but since the purpose here is 
to illustrate eigenvalues as functions of a parameter, only a brief description will be given. 


Modes are numbered by the real parts of their eigenvalues. This numbering is con- 
sistent with the number of beats or changes of 7 in the phase of the depth functions. Thus 
the eigenvalue of a mode numbered | ina profile consisting of only the upper duct lies 
exactly over the eigenvalues of a mode in the double duct in figures 8 and 9, but this mode 
in the double duct changes number each time it crosses the real part of another mode. The 
depth function actually gains an additional beat each time this happens. The background of 
modes that are being crossed consists of the higher order, untrapped modes associated with the 
lower duct. 


Mode 2, of the upper duct only, does not have a single mode in the double duct that 
overlies it exactly. Instead, a mode attempts to follow it at frequencies above 1350 Hz. 
Below this frequency, successive modes follow its path for short intervals. This interplay 
between modes occurs when mode 2 of the upper duct is in some sense equally as untrapped 
as the modes associated with the lower duct. 


The imaginary parts of the modes follow similar patterns; but because the mode 
numbering is not determined by the imaginary parts, the mode numbers sometime jump 
from one line to another. An important feature of these two plots is that if the real parts of 
the eigenvalues cross, the imaginary parts do not; and vice versa. Thus two eigenvalues do 
not tend to become equal at a point which would make them degenerate. 


The mode follower program will tend to follow the continuous curves. Thus if start- 
ed in the right direction on mode 59 at 1450 Hz, it will follow along the continuous mode 
which becomes successively mode 58, 57, 56, and 55. 


Figures such as 8 and 9 can be drawn by computing the eigenvalues at a sufficient 
number of frequencies to determine the lines. The mode follower does this for a given 
eigenvalue while adjusting the step size so the mode will not be lost, or so the program will 
correctly follow the mode. The step size is permitted to become large where the eigenvalue 
can be approximated by a parabolic curve, but it shortens when extrapolation to the next 
point becomes less accurate. 
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Figure 7. A five-layer approximation to the sound speed 
of a surface duct overlying a refractive duct. 
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Figure 8. The real part of the eigenvalues as continuous 
functions of frequency for the profile of figure 7. 

Some mode numbers are given. The first three modes 
for the surface duct only (SD) are shown as broken 
lines. That for mode 1 coincides with an existing line. 


37) 


MODE ATTENUATION COEFFICIENT, dB/km 


1150 1200 1250 1300 1350 1400 1450 


FREQUENCY, Hz 


Figure 9. The imaginary parts of the modes whose 
real parts are shown in figure 8 expressed as mode 
attenuation. To avoid confusion, they are not 
shown across the full frequency interval. 


When frequency is the variable parameter, the group velocity of the mode can be 
computed easily since a numerical derivative can be computed. Group velocity is given by 
the relationship 


Cy = dw/d (Re k) 


Ile 


Aw/A (Re k) 


IR 


-Afv2/f Av , 


where k is the horizontal wave number of the mode and v is the real phase velocity. The 
mode follower prints this value out at each step, along with the eigenvalue. 


IMPLEMENTATION OF THE MODE FOLLOWER 


The mode follower was originally implemented for a two-layer normal mode which 
differed from the n-layer program in that the derivative dG/dv of the characteristic equation 
was evaluated along with G. The iteration for roots of G was thus Newton-Raphson and is 
given by the relationship 


Vi+1 = Vi = G/G' " (86) 
This is simpler than the secant iteration of eq (15), in which G’ must be evaluated numeri- 


cally. Because of the simpler iteration, an effective scheme for mode following was available. 
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Since considerable effort was required to develop a similar scheme for the n-layer case, the 
two-layer mode follower will be briefly described to serve as an introduction to the n-layer 
case. 


The two-layer mode follower employs one iterative step of eq (86) at each point 
where G is evaluated. Thus, a root that is inexact but sufficiently exact is obtained. The 
original estimate is obtained by extrapolating from the three most recent roots. If this esti- 
mate is sufficiently close to the true root, the single iterative step will make a small correc- 
tion, G/G’, that will bring the estimate very close to the true root. By using the size of this 
correction to control the step size, the program is self-regulating. The program works well 
when a permissible value of G/G’ of 10-6 to 10-4 is used. Outside this interval the step size 
is either doubled or halved. 


The multiple-layer program differs from this in several details. The extrapolation 
from the previous three points is done not only for the phase velocity but also for the nu- 
merical derivative, 


D7! = Av/AG. 
Lagrange three-point interpolation is used, given by the form 


_ V(Ky) & = X59) (K = XZ) W(Kg) (K- x) K- XZ) VOR) (KX - x1) &- XQ) 
CTS )ICoES 7) MCS ES7)IC7 ES 7) MMCSEES UICT EST) 
(87) 


where x is the new value of the parameter that is being varied (usually frequency) and x], 
x7, and x3 are the three previous values, x; being the most recent. To extrapolate the 
derivative, v is replaced by D-! in eq (87). Both quantities are complex numbers. 


The determinant is now evaluated at this new phase velocity to give a value Gg. 
Next a corrected value of phase velocity, vg, is obtained: 


Yo = ¥-Gg Ds!” (88) 


In the two-layer case, the size of the correction, Gp"! , is used to control the step size. 
Because the numerical derivative is less precise, we evaluate G once more at this new posi- 
tion, obtaining G;. A new numerical derivative is next calculated: 


Do! = (o-v)/(G, -Gp) . 


This derivative is now compared with the extrapolated value to determine whether the step 
size should be changed. To do this an error 


E = |1-Do/DI? 


is computed. Good results have been obtained by keeping E between 1075 and 1072. IfE 
becomes larger than this, the step size is halved and the extrapolation is tried again. Should 
halving the step size five times fail to obtain a value of E less than 10-2, the mode is pre- 
sumed to be lost and the program halts. 
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If E is less than 10-2, the step is successful and the stored values are updated for the 
next step. Before v is stored, though, the iterative step of eq (88) is applied one more time 
to obtain a more precise value of v. This requires little extra effort because the quantities 
G and Do7! are already available. 


If the error E is less than 1079, the next step size is doubled. 


It is possible for the extrapolation to be too successful. That is, if v is very near the 
true root, Gg and Av will be very small. The numerical derivative may then be inaccurate. 
Therefore, before the error term E is computed, a quantity 


F = iv/Avi2 


is computed. If F is greater than 1028 the extrapolated derivative is used rather than the 
computed derivative and the program proceeds to the next step. If F is greater than 1034, 
the step size is doubled before proceeding to the next step. 


The other principal part of the program is the initialization which must evaluate v at 
three values of x to obtain the numbers needed for the first extrapolation, eq (87). 


INPUT AND OUTPUT 


The first input card contains the maximum number of steps allowed, the limits 
applied to E and F, and keys which control both the amount of detail in the printout and 
whether the profile parameters are to be read in or retained from the previous run. Default 
values are supplied when these items are left blank. Next the profile parameters are read in. 
These are an older style and only permit specification of the absorption loss at the top of a 
layer. The sound speed gradient is assumed to be real at the top of any layer. 


A final card indicates which variable — frequency, sound speed, depth, gradient, 
absorption, or density — will be varied, by specifying a number called nx in the program, 
from 1 to 6. The next number, ny, specifies which layer the variable will be in. This layer 
number is not needed if frequency is selected. A third number, nz, indicates, if zero, that 
the profile will remain continuous as the selected parameter is varied. If nz is not zero, the 
selected parameter moves alone without a compensating motion in other profile parameters. 
The card next gives the initial and final value of the parameter to be varied and the initial 
step size. Finally, the particular mode to be followed is indicated by giving an approximate 
phase velocity and an initial step size. These must be chosen such that the subsequent itera- 
tion will converge on the correct mode. 


The principal output of this program is the print statement at line 314. Each line of 
output contains the value of the parameter being varied, the complex phase velocity, the 
determinant G, the derivative Dale the error term E, the mode attenuation, the mode group 
velocity (if frequency is the parameter being varied), and the step number. After the final 
step, the profile in its final form is printed out. 
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CONCLUSIONS 


1. An effective program for computing propagation loss in a layered ocean by nor- 
mal modes has been developed. Complete documentation for the program is contained 
herein. 

2. Sediment layers are modeled as fluids in which densities, sound speeds, and ab- 
sorption can be specified. This permits a complete wave solution for bottom reflected 
sound energy. 

3. A continued fraction technique for evaluating asymptotic series is shown to give 
superior results in evaluating the auxiliary functions required in this program, the modified 
Hankel functions of order 1/3. 


4. A mode follower program given here is useful in tracing eigenvalues. Such traces 
are needed to understand the eigenvalue structure. 


RECOMMENDATIONS 


1. Improve the mode locating ability of this normal-mode program to make it self- 
contained. It currently requires user interaction to locate eigenvalues. 


2. Investigate methods to incorporate the effect of rough boundaries into this 
program. 
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APPENDIX A: NORMAL MODE PROGRAM IN FORTRAN 


This program consists of the main program and seven subroutines. The main pro- 
gram handles the input and output and performs much of the computation. This includes 
profile preparation, mode search, determination of depth function coefficients, normaliza- 
tion, computation of depth functions, and summation of modes. Auxiliary functions are 
performed by the subroutines SETUP and DET, which set up the determinant, then evaluate 
it. This is the determinant from which eigenvalues are determined. The subroutine HZERO 
determines the Hankel functions of order zero, second type, which gives the range depend- 
ence of the modes. Only a single term of the asymptotic expansion is needed for this 
function. 


Subroutine HANKEL evaluates the modified Hankel functions of order 1/3, by 
which the depth dependence of the modes is expressed. The majority of computing time is 
usually expended in this subroutine. Subroutine CFR is used by subroutine HANKEL to 
evaluate continued fractions. Subroutine RCOEF evaluates and prints reflection coeffi- 
cients when they are requested. 
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OODNOuUDBWNH— 


Cc 


THIS IS THE MAIN PART OF NLAYNM 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION LAMBDA, LAMBDI 

INTEGER COL 

REAL R ATTEN, T RE, RX 

DIMENSION LOSPCH(5,25) 

COMMON /HAN/ H2R,H21,H1R,H11,H2PR,H2PI,H1PR,HI1PI,R 

COMMON/INPUT/ Z(12), N, OMEGA, V, VI, CON(12), GSQ(12), 

1 CAY(12), LAMBDA, LAMBDI, G(12) 
2,RHO(12), GI(12), G SQI(12), CAYI(12) 

COMMON /EXPO/ EXSUM, CNTR, RATIO(25) 

COMMON/DETMNT/ A(25,4), Q(25,4) 

COMMON/PARTS/ ZT(12),ZT1(12),ZB(12),ZB1(12) 

COMMON/REFL/ AF(12,200), AG(12,200), BF(12,200), BG(12,200), 
2 EIGEN(350), EIGENI(350), B (25,4), BI(25,4), CB(12), CBI(12), 
3 CAYSQ(12), CAYSQI(12), NN 

DIMENSION D(350), DI(350), F(100), F1(100), HZERO2(350), 

* DA(350), SRES(350), GAMMAI(12), BLPK(12), 
4 HZER21(350), DPK(12), GCU(12), GCUI(12), C1(12), 
3 PHASE v(350), PHASI V(350), UU(2000), UUI(2000) 

COMMON /LIMIT/ TLIM, EXPONT, SLIM 

DIMENSION LOSS(101) 

DIMENSION C(12), DEPTH(52), DBLOSS(350), COL(120), 
1CONTR(10), EF(2), FMAG(350), FANG(100), 
2GAMMA(12),JSMBL(10) ,JCOUNT(5),JUCOU(5),LEVEL(41),PLEV(5) ,RLOSS(100) 


3 , RLOS(101),RECVRS(51),TEST(3), ING(11) 
EQUIVALENCE (FF,EF(1)),(DEPTH(1),SOURCE), (DEPTH(2),RECVRS(1)), 
1 (RLOS(2),RLOSS(1)) 


COMMON /AHZERO/ HZEROR,HZEROI 

DATA ( CONTR(I), 1=1,4) /110.00,95.D0,80.D0,-1000.D0/, 

1 (JU SMBL(I), I=1,3) /1H1,1H*,1H8/, 

*(ING(I), I = 1,10)/1HO,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/, 
2 (TEST(I), I1=1,3) /.2D0,1.D0,5.D0/ 

TLIM = 25. 

SLIM = -8.54988 


READ IN PARAMETERS 


1 


READ 11, K1, K2, K3, K4, K5, K6, K7, K8, KY 


C KEYS& 1-DEPT FN PRINTOUT, 2-LOSS PRINTOUT, 3-REFLECTION COEFF PRINTOUT 


Cc 


4- 


11 


13 


434 


CHANGE CONTOURS, 5-CONTINUE MODES 
FORMAT (1014) 
FRINT 13, h1, KZ, K3, K4, KS, K6, K7, K8, hKY 
FORMAT (GH1IKEYS , 1014) 
i? (Ci ‘oka 7) Cok ves 
READ 11, M 
READ 434, (SRES(I), I = 1,M) 
FORMAT (5D16.7) 
IF (K8 .NE. 1) GO TO 8 
READ 20, T LIM 
SLIM = -DCBRT(TLIM**2) 
PRINT 30, TLIM, SLIM 
EXPONT = DEXP((TLIM + TLIM) / 3.) 
K6 = K6 + 1 
INSUR = 0 
IF (K2 .LT. 10) GO TO 16 
K2 = K2 - 10 
INSUR = 1 
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16 


17 


20 


18 
19 
14 


30 


MPCH = 0 

1F (K5 .LT. 10) GO TO 17 

K5 = K5 - 10 

MPCH = 1 

IF (K4 .NE. 1) GO TO 3 

READ 20, (CONTR(I), I = 1,9) 
CONTR(10) = -1000.D0 

READ 4, (JU SMBL(I), I = 1,9) 
FORMAT (9A1) 

READ 10, N, FREQ 

FORMAT (12, D10.1) 

IF (N.EQ.0) GO TO 999 

PRINT 12, N, FREQ 

FORMAT (13, 8H LAYERS, ,F10.1, 
READ 20, (Z(I), I=1,N) 


PRINT30, (Z(I), I=1,N) 
READ 20, (C(I), I = 1,N) 
FORMAT (8D10.4) 

PRINT3O, (C(I), I = 1,N) 
READ 20, (CB(I), I = 1,N) 
PRINT30, (CB(I), I = 1,N) 
READ 20, (GAMMA(T), I = 1,N) 
PRINT30, (GAMMA(I), I = 1,N) 
READ 20, (DPK(I), I = 1,N) 
PRINT 30, (DPK(I), I = 1,N) 
READ 20, (BLPK(1), I = 1,N) 
PRINT3O, (BLPK(I), I = 1,N) 
READ 20, (RHO(I), I = 1,N) 
PRINT 30, (RHO(I), I = 1,N) 
IF (FREQ .GT. 0.) GO TO 18 
FREQ = - FREQ 

ATTEN = 0. 

GO TO 19 


F SQ = (FREQ / 1000.)**2 

ATTEN -1 * F SQ / (1. + F SQ 
ATTEN = ATTEN * 1.0936 

PRINT 14,ATTEN 

FORMAT (8H ATTEN = ,G10.5, 5HD 
ATTEN = ATTEN / 1000.D0 

FORMAT (9F14.5) 


C COMPLETE PROFILE 


31 


32 


34 


36 


DO 33 I = 1,N 

IF (RHO(I) .EQ. 0.) RHO(I) = 1 
IF (CB(1) .NE. 0.) GO! TO) 31 
CB(I) = C(I+1) 

IF (C(1) .NE. 0.) GO TO 32 
C(I) = CB(I-1) 

IF (DPK(I) .GE. 0.) GO TO 34 
GI(I) = CBI(I-1) 

GO TO 36 

CI(1) = 0. 

IF (DPK(I) .EQ. 0.) GO TO 36 
T = 27287.52708 / DPK(1) 

CI(1) = T - SQRT((T - C(1)) * 
CBI(I) = 0. 

IF (GAMMA(I) -NE. 0.) GO TO 38 
**BOTH SOUND SPEEDS GIVEN 


3H HZ) 


) + 40. 


B/KM ) 


02 


(T + C(1))) 
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* F SQ / (4100. + F SQ) 


114 IF (BLPK(I) .LE. 0.) GO TO 37 


115 1 = 2/28/.92/0%8 / BLPK(I) 

116 CBI(I) = T - SQRT ((T - CB(I)) * (T + CB(I))) 

117 37 T = C(I) * (C(1)*#*2 - 3. * CI(1)**2) 

118 TI = CI(I) * (3. * C(1)**2 — CI(1)**2) 

119 IF (BLPK(I) .LT. 0.) GO TO 39 

120 TEMP = CB(I)**2 - CBI(I)**2 

121 TEMP1 = 2. * CB(I) * CBI(1) 

122 DENOM = TEMP**2 + TEMP1**2 

123 TEMP = TEMP / DENOM 

124 TEMP1 = -TEMP1 / DENOM 

125 GAMMA(I) = 0.5 * (C(I) -— (T * TEMP - TI * TEMP1)) / 
126 * (Z(I+1) - Z(1)) 

127 GAMMAI(I) = 0.5 * (CI(I) - (T * TEMP1 + TI * TEMP)) / 
128 * (Z(I+1) - Z(1)) 

129 IF (I .EQ. N) GO TO 27 

130 GO TO 33 

131 c **SPECIAL CASE, GRADIENT REAL NUMBER 

132 39 IF (CI(1I) .EQ. 0.) GO TO 42 

133 TEMP = CB(I1)**2 

134 TEM = TEMP**2 

135 TEMP1 = CI(1) 

136 COEF1 = CI(1) 

137 COEF2 = 2. * TEMP * CI(I) + TI 

138 COEF3 = 2. * T * CB(I) 

139 COEF4 = TEM * CI(1I) - TEMP * TI 

140 OLDFN = 1.020 

141 DO 41 Jv = 1,10 

142 FN= (((COEF1 * TEMP1) + COEF2) * TEMP1 + COEF3) * TEMP1 + COEF4 
143 FP = ((4. * COEF1 * TEMP1) + 2. * COEF2) * TEMP1 + COEF3 
144 TEMP1 = TEMP1 - FN / FP 

145 IF (FN .GE. OLDFN) GO TO 43 

144 OLDFN = FN 

147 41 CONTINUE 

148 43 CBI(I) = TEMP1 

149 GAMMA(I) = .5 *(.5 * (CI(I) *(TEMP - CBI(I)**2) = TI) / 
150 * (CB(I) * TEMP1) + C(I)) / (Z(1+1) - Z(1)) 

151 GO TO 28 

152 42 GAMMA(1) = C(I) - T / (CB(1I)**2 * (Z(I+1) - Z(I)) * 2.) 
153 GO TO 33 

154 C **SOUND SPEED AND GRADIENT GIVEN 

155 38 IF (I .£Q. N) GO TO 33 

156 T = C(I) * (C(I)**2 - 3. * CI(1)**2) 

157 TI = CI(1I) * (3. * C(1)**2 - CI(1)**2) 

158 IF (BLPK(1) .EQ. 0.) GO TO 29 

159 IF (BLPK(1) .LT. 0.) GO TO 28 

160 TEMP = (BLPK(I) / 54575.05416)**2 / (Z(I+1) - Z(1)) * 0.5 
161 T = T * TEMP 

162 TI = TI * TEMP 

163 TEMP = .5 * C(I) / (Z(I+1) - Z(1)) - GAMMA(I) + T 

164 T = -(TI - SQRT( TI * TI + T * TEMP)) / T 

165 CB(I) = 54575.05416 * T / BLPK(I) / (1. + T * T) 

166 CBI(I) = CB(I) / T 

167 GO TO 37 

168 Cc **SPECIAL CASE, GRADIENT REAL NUMBER 

169 28 TEMP = C(I) -— 2. * GAMMA(I) * (Z(1+1) - Z(1)) 

170 TEM = TEMP**2 + CI(1)**2 


46 


XRE = (T * TEMP + TI * CI(1)) / TEM 
XIM = -(T * CI(I) - TI * TEMP) / TEM 
TEM = XRE**2 + XIM**2 


CB(I) = SQRT((XRE + SQRT(TEM)) * .5) 
CBI(I) = .5 * XIM / CB(I) 
GAMMAI(1) = 0. 
GO TO 33 
29 TEMP = C(I) - 2. * GAMMA(I) * (Z(I+1) - Z(I)) 
CB(I) = SQRT (T / TEMP) 
GAMMAI(I) = .5 * (CI(I) - TI / CB(1)**2) / (Z(I+1) - 2(1)) 
GO TO 33 
27 N=N-1 
33 CONTINUE 
C COMPUTE USEFULL QUANTITIES 
PRINT 58 
58 FORMAT (7X,6H RE M ,8X,6H IMM ,9X,5H L/KM,8X,6H RE C ,8X, 
* 6H IM C ,5X,12H RE C BOTTOM,4X,12H IM C BOTTOM,10X,9H GRADIENT ) 
OMEGA = 6.283185307D0 * FREQ 
DO 40 I = 1,N 
TEMP = C(1)**2 + CI(I)¥**2 
CAY(1) = OMEGA = C(I) / TEMP 
CAYI(1) = -OMEGA * CI(I) / TEMP 
CAY SQ(I) = CAY(I)**2 - CAYI(1)**2 
I 


CAY SQI(I) = 2.D0 * CAY(I1) * CAYI(I) 

TEMDR = -2. * (GAMMA(I) * CAY SQ(I) - GAMMAI(I) * CAY SQI(I)) 
TEMDI = -2. * (GAMMA(I) * CAY SQI(I) + GAMMAI(I) * CAY SQ(T)) 

G CU(I) = (TEMDR * C(I) + TEMDI * CI(I)) / TEMP 

G CUI(I1) = (TEMDI * C(I) - TEMDR * CI(I)) / TEMP 

TEM1 = DCBRT(-DSQRT( GAMMA(I)**2 + GAMMAI(1)**2) * 2.*OMEGA**2) 


TEM1I = DATAN (ABS(GAMMAI(1) / GAMMA(I)))/ 3. 
CRTG = TEM1 * DCOS(TEMI1T) 

CRTGI = TEM1 * DSIN(TEM11) 

IF (GAMMA(I) .LT. 0.) CRTG = -CRTG 

IF (GAMMAI(I).LT. 0.) CRTIGI = -CRTGI 

G(I) = (C(I) * CRTG + CI(1) * CRTGI) / TEMP 


GI(I) = (C(1) * CRTGI - CI(I) * CRTG) / TEMP 
CON(I) = G(1) * C(I) - GI(1) * CI(1) 

CON(I) = OMEGA**2 / CON(1)**2 

XMI = -GI(I) * (Z(1+1) - Z(1)) 


XM = -G(I) * (2Z(I+1) - Z(1)) 
DPK(I) = -8686.D0 * CAYI(1) 
PRINT 30, XM, XMI, OPK(I), C(I), CI(I), CB(I), CBI(1I) 
* ,GAMMA(I), GAMMAI (TIT) 
G SQI(1) = 2. * G(I) * GI(I) 
40 G SQ(I) = G(I)**2 - GI(I)**2 
C FIND MODES 
NXTRA=0 
IU FLAG=0 
NN = NN + 1 
IF (KS .EQ. 1) GO TO 15 
DO 50 NN = 1,350 
15 IF (IU FLAG .EQ. 1) GO TO 53 
52 IF (NXTRA .GT. 0) GO TO 44 
READ 60, V,VI,STEP,STEPI,NXTRA 
60 FORMAT (4D10.4,110) 
IF (NXTRA .GE.0) GO TO 62 
V=V+VI * 1.0710 


47 


62 

14 
Cc 

14 


55 


54 


56 


45 


47 


44 


70 


80 


2 


VI = STEP + STEPI * 1.D-10 
GU fu 85 

IF (V) 142,301,70 

IF(STEP) 44,444,143 


SEARCH FOR MODE 


3 


46 


53 


SUZESm— ei 

SIZE2=0 

IJU FLAG=1 

V=-V 

IF (NXTRA) 55,55,54 

NXTRA = 20 

XTRA = NXTRA 

HOP = (STEP - V) / XTRA 

HOPI=0. 

IF(STEPI.NE.O.) HOPI=(STEPI-VI)/XTRA 
DO 47 IU = 1,NXTRA 

CALL SETUP 

DET = VEL 

DETI = VELI 

CALL DETNT(N,VEL,VELI) 

DELTA = VEL 

MEL TT - Weit 

SIZE = DELTA*DELTA + DELTI*DELTI 
BIRUNIT (SiGe envi uViLi Syli Zier ue VIEL AVE ST: 
FORMAT (2F12.3, 3017.5) 

IF ((SIZE2.LT.SIZE3).AND.(SIZE.GT.SIZE2)) GO TO 45 
SIZE3=SIZE2 

SIZE2=SIZE 

V = V + HOP 

VI=VI+HOPI 

GO TO 47 

Ve= Ve = = HOP 

TEMP = HOP / (SIZE - SIZE2) 

DELTI = TEMP * (DET * VELI - DETI * VEL) 


TEMP = .5D0 * (SIZE3 - SIZE) / (SIZE3 + SIZE - SIZE2 - SIZE2) 


DELTA = HOP * TEMP 

IF(HOPI.EQ.0) GO TO 49 

VI=VI-HOPI 

DELTAI=HOPI*TEMP 

GO To 49 

CONTINUE 

IU FLAG=0 

NXTRA=0 

GO TO 52 

SIZE 2=-1. 

SIZE=0 

GO TO 46 

NXTRA = NXTRA - 1 

V = 3. * (PHASE V(NN-1) - PHASE V(NN-2)) + PHASE V(NN-3) 
VI = 3.* (PHASI V(NN-1) - PHASI V(NN-2)) + PHASI V(NN-3) 
STEP = (PHASE V(NN-1) — PHASE-V(NN-2)) * .0001 
CALL SETUP 

CALL DETNT(N,DET,DETI) 

FORMAT (/, 2020.11, 4013.4) 

VEL = DET 

VELI = DETI 

DELTA = STEP 


48 


285 DELTI = STEPI 


286 1F (DELIA .NE.O.) GO TO 49 

287 IF (DELTI .EQ.0.) DELTA = .01 

288 49 SIZE2 = 100. 

289 RX = DET**2 + DETI**2 

290 IF (K6 .LT. 3) PRINT 80, V, VI, DET, DETI, SIZE, CNTR 
291 vu =0 

292 oH dads 

293 IF (Uv .GT. 15) GO TO 51 

294 V = V + DELTA 

295 VI = VI + DELTI 

296 IF (VI) 83,84,85 

297 83  ODELTI = DELTI - VI 

298 84 VI = 1.D-18 

299 85 CALL SETUP 

300 NNN = N+N—- 1 

301 DO 82 IA = 1,NNN 

302 DO 82 1B = 1,4 

303 BI(IA,1B) = Q(IA,1B) 

304 82 B(IA,IB) = A(IA,1B) 

305 CALL DETNT(N,DET,DETI) 

306 IF (K6 .NE. 1) GO TO 72 

307 71. PRINT 81, V, VI, DET, DETI, SIZE, CNTR 
308 81 FORMAT (2020.11, 4013.4) 

309 72 IF (NXTRA .LT. 0) GO TO 51 

310 TEMNR = DET * DELTA —- DETI * DELTI 
311 TEMNI = DETI * DELTA + DET * DELTI 
312 TEMDR = VEL - DET 

313 TEMDI = VELI - DETI 

314 TEMDEN = TEMDR*TEMDR + TEMDI*TEMDI 
315 IF (TEMDEN .EQ. 0.) GO TO 51 

316 TEMRNU = TEMNR*TEMDR + TEMNI*TEMDI 
317 TEMINU = TEMNI*TEMDR - TEMNR*TEMDI 
318 DELTA = TEMRNU/TEMDEN 

319 DELTI = TEMINU/TEMDEN 

320 C * * * THE NEXT CONSTANT DEPENDS ON WORD LENGTH AND SIZE OF PHASE VELOCITY * * 
321 IF (ABS(DELTA) .LT. 1.0-14) GO TO 51 
322 SIZE = DELTA*DELTA + DELTI*DELTI 

323 IF ((SIZE.GT.SIZE2).AND.(U.GT.3)) GO TO 51 
324 G2F ST ZEQNamS ize non 

225 VEL = DET 

326 VELI = DETI 

327 GO TO 48 

328 C FIND DEPTH FUNCTIONS 

329 51 IF (INSUR .EQ. 0) GO TO 61 

330 TRE = (DET**2 + DETI**#2) / RX 

331 IF (TRE .LT. 1€-10) GO TO 61 

332 PRINT 998, NN, TRE 

333 998 FORMAT (5H MODE ,14,23H FAILED TO CONVERGE -- , E9.2) 
334 GO TO 999 

335 61 IF (MPCH .EQ. 0) GO TO 63 

336 IF (NXTRA .LT. 0) GO TO 63 

337 TEM1 = V * 1.04 

338 COL(1) = TEM1 

339 TEMP = COL(1) 

340 COL(2) = (TEM1 - TEMP) * 1.D10 

341 TEM1 = VI * 1.04 


49 


64 
63 


74 
73 


96 


97 


1 
98 


95 


Cc 


1 


1 
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COL(3) = TEM1 

TEMP = COL(3) 

COL(4) = (TEM1 - TEMP) * 1.010 
COL(5) = -NN 

PUNCH 64, (COL(I), I = 1,5) 
FORMAT (5110) 


AF(1,NN) = A(1,3) 
AG(1,NN) = Q(1,3) 
BF(1,NN) =-A(1,4) 
BG(1,NN) =-Q(1,4) 
PHASE V(NN) = V 

PHASI V(NN) = VI 


IF (K6 .EQ. 1) GO TO 73 

PRINT 81, V, VI, DET, DETI, SIZE, CNTR 

(= SING iat 

IF (LL-1) 95,96,97 

i BO 

GO TO 98 

DO 110 Jv = 2,LL 

Woo Wd od 2.2 

TEMNe A(1,2)*AF( U-1 ,NN) - Q(1,2)*AG( JU-1 ,NN) 

TEMNI = Q(1,2)*AF( U-1 ,NN) + A(I,2)*AG( JU-1 ,NN) 

TEMDR = A(I,3)*A( I+1 ,4) - Q(1,3)*Q( [+1 ,4) - 
A(I,4)*A( I+1 ,3) + Q(1,4)*Q( I+1 ,3) 

TEMDI = Q(1,3)*A( I+1 ,4) + A(1I,3)*Q( I+1 ,4) =—- 
Q(1,4)*A( 141 ,3) - A(1,4)*Q( I+1 ,3) 


TEMDEN = TEMDR*TEMDR + TEMDI*TEMDI 
TEMRNU = TEMNR*TEMDR + TEMNI*TEMDI 
TEMINU = TEMNI*TEMDR —- TEMNR*TEMDI 


TEMP = TEMRNU / TEMDEN 
TEMPI = TEMINU / TEMDEN 


BF(JUJ,NN) = -(TEMP+A( I+1 ,4) - TEMPI*Q( I+1 ,4)) 
BG(JUJ,NN) = -(TEMPI*A( I+1 ,4) + TEMP*Q( I+1 ,4)) 
AG(JU,NN) = TEMPI*A( I+1 ,3) + TEMP*Q( I+1 ,3) 

AF(JUJ,NN) = TEMP*A( I+1 ,3) - TEMPI*Q( I+1 ,3) 

TEMNR = -(A(1I+2,2) * AF(LL,NN) -Q(1+2,2) * AG(LL,NN) ) 
TEMNI = -(Q( I+2 ,2)*AF(LL,NN) + A( I+2 ,2)*AG(LL,NN) ) 
TEMDEN = A( I+2 ,3)*A( I+2 ,3) + Q( I+2 ,3)*Q( I+2 ,3) 
TEMRNU = TEMNR*A( I+2 ,3) + TEMNI¥Q( I+2 ,3) 

TEMINU = TEMNI*A( I+2 ,3) - TEMNR*Q( I+2 ,3) 

BF(N,NN) = TEMRNU / TEMDEN 

BG(N,NN) = TEMINU / TEMDEN 

AF(N,NN) = 0. 


AG(N,NN) 


= 0. 
FIND NORMALIZING FACTOR 


1 


D(NN) = 2.12429296D0 * RHO(1)**3 / G(1) 


DI(NN) = 0. 

DO 111 I = 2,N 

TEMRSP = AF( I-1 ,NN)*B( 2*I-2 ,2) - AG( I-1 ,NN)*BI( 2*I-2 ,2) 
BF( I-1 ,NN)*B( 2*I-2 ,1) - BG( I-14 ,NN)*BI( 2*I-2 ,1) 


TEMISP = AG( I-1 ,NN)*B( 2*I-2 ,2) + AF( I-1 ,NN)*BI( 2*I-2 ,2) 
BG( I-1 ,NN)*B8( 2*I-2 ,1) + BF( I-1 ,NN)*BI( 2*1-2 ,1) 

AX1 = TEMRSP*TEMRSP - TEMISP*TEMISP. 

AX1I = TEMRSP * TEMISP 

AX11 = AX1I + AX1I 

(G(1-1)**2 + GI(I-1)**2) 

G(I)**2 + GI(1)**2 


4 
m 

= 

oO 

Po) 
wou 
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Cc 


11 


TEMP = (RHO(I-1) 
TEM1 =(ZB (1-1) 
* -(ZT (1) * G(T) 
TEM1I =(ZBI(I-1) 
* -(ZTI(1) * G(1) 


RHO(I)) / TEMDI 


ZTI(1) * GI(1)) * TEMP 


> 4p TOSS 


ZT (1) * GI(1)) * TEMP 


G(I-1) + ZBI(I-1) * GI(I-1)) / TEMDR 
G(I-1) - ZB (I-1) * GI(I-1)) / TEMOR 


»NN)*BI( 2*I-1 


TEMRSP = AF( I-14 ,NN)*B( 2*I-1 ,2) - AG( I-71 
1 BF( I-1 ,NN)*B( 2*I-1 ,1) - BG( I-1 ,NN)*BI( 2*I-1 
TEMISP = AG( I-1 ,NN)*B( 2*I-1 ,2) + AF( I-1 


»NN) *BI( 2*I-1 


1 BG( I-1 ,NN)*B( 2*I-1 ,1) + BF( I-1 ,NN)*BI( 2*I-1 


AX2 = TEMRSP*TEMRSP - TEMISP*TEMISP 
AX2I = TEMRSP * TEMISP 
AX2I = AX2I + AX2I 


TEMDR = RHO(I-1) / (G CU(I-1)**2 + G CUI(I-1)**2) 


TEMDI 


RHO(I ) / (G CU(I )*¥*2 + G CUI(I 
TEM2 = G CU(I-1) * TEMOR - G CU(I) * TEMOI 


)**2) 


TEM2I = G CUI(I) * TEMDI - G CUI(I-1) * TEMDR 
TEMR1 = AX1*TEM1 - AX1I*TEMII 

TEMI1 = AX11*TEM1 + AX1*TEMII 

TEMR2 = AX2 * TEM2 - AX2I * TEM2I 

TEMI2 = AX2I * TEM2 + AX2 * TEM2I 

D(NN) = D(NN) + TEMR1 / RHO(I-1) + TEMR2 


DI(NN) = MTONN) + TEMI1 / RHO(I-1) + TEMI2 


1 CONTINUE 


IF ( K1 .GT. 3) DA(NN) = DSQRT((D(NN)**2 + DI(NN)**2) * FREQ / 


* PHASE V(NN) ) 
EIGEN(NN) = LAMBDA 
EIGENI(NN) = LAMBDI 
IF (K6 .GT. 2) GO TO 131 


oO 

wou 
= 
Lee =s 
— 


= 1,N 


+ 


1 
) = SNGL(ZT(I)) * 100. 
1 
) = SNGL(ZTI(1)) * 1000. 
1 
) = SNGL(ZB(I)) * 100. 
1 


(=) 


+ 


oO oO 
oreo w on 
=~ 


=~ 
RRAOr OT KKANO 
+ 


RFOROFOCFVOKEK 


+ 


COL(K) = SNGL(ZBI(I)) * 1000. 


112 CONTINUE 


PRINT 130, (COL(I), I=1,L) 
PRINT 130, (COL(I), I1=25,kK) 
130 FORMAT (4H Z = , 11(16,15)) 
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30 


M=N+N 
PRINT 132, (RATIO(I), I = 1,M) 
2 FORMAT (11(1X,2F5.3)) 
1 OB LOSS(NN) = - LAMBDI * 8686.D0 


PHINV = V * PHASE V(NN-1) /((V - PHASE V(NN-1))* FREQ) 
PRINT 109, NN,EIGEN(NN),EIGENI(NN),O(NN),DI(NN) ,PHINV,OB LOSS(NN) 
9 FORMAT (3H N=,I5,10H LAMBDA =,2E15.7,4H D= 


* 12H INT RANGE = ,F8.0, 6H L/K =, F8.5) 
IF (K3 .EQ. 0) GO TO 50 
CALL RCOEF (K3) 
CONTINUE 
READ IN SOURCE AND RECEIVERS DEPTHS 
1 NRT = NR 
NR = 0 


Sul 


»2E15.7, 


NN = NN - 1 


K1P1 = K1 + 1 

IF (Ki .NE. 3) GO TO 321 
NR = NRT 

GO TO 501 


321 READ 20, SOURCE 
320 NR = NR + 1 
READ 20, RECVRS(NR), FINAL, STEPP 
IF (NR.GT.50) GO TO 300 
350 IF (RECVRS(NR) .EQ.0.) GO TO 300 
310 IF (FINAL .EQ.0.) GO TO 320 
330 RECVRS(NR+1) = RECVRS(NR) + STEPP 
IF (RECVRS(NR+1) .GT. FINAL) GO TO 320 
340 NR = NR + 1 
IF(NR .GT. 50) GO TO 300 
GO TO 330 
300 PRINT 303 
303 FORMAT (/21H SOURCE AND RECEIVERS ) 
PRINT 21,(DEPTH(I), I = 1,NR) 
21 FORMAT (8F10.2) 
C COMPUTE DEPTH FUNCTIONS 
DO 500 I = 1,NN 
Loc = 1 
DO 305 J = 1,NR 
IF ((JU .EQ. 1) -AND. (K1 .GT. -5)) GO TO 305 
LCTR = 0 
380 IF((DEPTH(U) .GE. Z(LOC)).AND.(DEPTH(JU) .LT. Z2(LOC+1)))GO TO 360 
371 IF (LOC .GE. N) GO TO 385 
370 LOC = LOC + 1 


GO TO 380 
385 IF (DEPTH(J) .GE. Z(LOC)) GO TO 360 
390 LOC=1 


LCTR=LCTR+1 
IF (LCTR .GT. 2) GO TO 305 


GO TO 380 

360 x1 = CAY (LOC) - EIGEN (1) 
X2 = CAY (LOC) + EIGEN (1) 
X3 = CAYI(LOC) - EIGENI(1) 
X4 = CAYI(LOC) + EIGENI(1) 


TEMP = X1 * X2 - X3 * X4 
TEMPT = X1 * X4 + X3 * X2 
TEMDEN = G SQ(LOC) **2 + G SQI(LOC)**2 
ZE = (TEMP * GSQ(LOC) + TEMPI * G SQI(LOC)) / TEMDEN 
ZEI = (TEMPI * GSQ(LOC) + TEMP * GSQI(LOC)) / TEMDEN 
TEM1 = ZE 
IF (ZE .GT. -7.5) GO TO 438 
= CAY(LOC) 
= CAYI(LOC) 
DO 437 K = 1,20 
TEMP = S¥*2 + T**2 
TEMPI = (EIGENI(I) * S - EIGEN(I) * T) / TEMP 
TEMP = (EIGEN(I) * S + EIGENI(I) * T) / TEMP 
ZE = ((1. + TEMP) * (1. - TEMP) + TEMPi**2) * CON(LOC) 
ZEI = -2. * TEMPI * TEMP * CON(LOC) 
ZR = ZE / -7.5 
IF (DABS(ZR-1.) .LT. 1.0-3) GO TO 438 
S = EIGEN(1) + (S - EIGEN(I)) / ZR 
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437 CONTINUE 
438 IF (G(LOC) .LT. 0.) GO TO 439 
ZE = G(LOC) * (DEPTH(JU) - Z(LOC)) + TEM1 
IF (ZE .GT. -7.5) GO TO 442 
F(J) = 1.0-12 
FI(J) = 0. 
GO TO 305 
439 ZE = G(LOC) * (DEPTH(J) - Z(LOC)) + ZE 
IF (ZE .GT. SLIM) GO TO 442 
F(J) = 1.D-12 
FI(J) = 0. 
GO TO 305 
442 ZEI =GI(LOC) * (DEPTH(JU) - Z(LOC)) + ZEI 
302 CALL HANKEL(ZE,ZEI,1) 


F(JUJ) =(AF(LOC,I)*#H1R - AG(LOC,1)*H1I + BF(LOC,I)*H2R - BG(LOC,1!) 


1 *H21) * RHO(LOC) 


FI(J) =(AG(LOC,I)+*H1R + AF(LOC,I)*H1I + BG(LOC,1)*H2R + BF(LOC,!I) 


1 *H21) * RHO(LOC) 
305 CONTINUE 
IF (K1 .—£Q. 2) GO TO 451 
fa) 3) G)si2 
451 PRINT 270, DEPTH(NR) 


270 FORMAT(7H1DEPTH ,F5.1,6X,3HE-8,17X,3HE-6,17X,3HE-4,17X, 3HE-2, 


* 17X,3HE O ) 
432 IF (K1 .LT. 4) GO TO 431 
IF (K1 .GT. 5) GO TO 433 
SRES(1) = (F(1)**2 + FI(1)**2) / DA(I) 
GO TO 500 
431 TEMDEN = D(1)*D(I) + DI(1I)*DI(1) 
TEMRE = F(1)*O(1) + FI(1)*DI(1) 
FD = TEMRE/TEMDEN 
FDI = (D(I) * FI(1) - DI(I) * F(1)) / TEMDEN 
433 DO 400 K = 2,NR 
vy ako y 
L = uJ * NN - NN# I 
IF (Ki .LT. 6) GO TO 435 
FF = SRES(1) * (F(K)**2 + FI(K)**2) / DA(T) 
GO TO 436 
435 FF = FD * F (K) - FDI * FI(K) 
FFI = FD * FI(K) + FDI * F(K) 
436 UU(L) = FF 
UUI(L) = FFI 
452 GO TO (400,410,420,400,400,400,400,400), KiP1 
C PLOT DEPTH FUNCTIONS 
420 DO 210 II = 1,120 
210 COL(II)= 1H 
DO 220 II= 20,100, 20 
220 COL(II)= 1HI 
FE = FF * FF + FFI * FFI 
IF ((FE.GT.1E-20).AND.(FE.LT.10000.)) GO TO 240 
GO TO 250 
240 INT = 100.D0 + 2.17147D0 * DLOG(FE) 
COL(INT) = 1H* 
GO TO 225 
250 COL(2) = 1H* 
225 PRINT 260, COL 
260 FORMAT (120A1) 
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c 


Cc 


GO TO 400 


PRINT DEPTH FUNCTIONS 
F MAG(J) = SQRT (FF * FF + FFI * FFI) 


410 
430 
440 
450 
170 
180 
400 
441 


500 


CALCULATE ATTENUATION AND READ IN RANGES 
IF ((K1 .EQ. 4).0OR.(K1 .EQ. 5)) PRINT 
(SRES(K), 


501 
551 


533 


562 


772 
561 


524 


563 


564 


560 
550 


* 


IF (FF) 430,440,450 


F ANG(J) = ATAN(FFI / FF) * 57.29577951D0 + 180.D0 


GO TO 400 

F ANG(J) = 90. 

GO TO 400 

F ANG(J) = ATAN(FFI / FF) 


FORMAT ( 10F12.4) 
FORMAT(/10E12.3) 
CONTINUE 

IF (K1.£Q.1) GO TO 441 
GO TO 500 

PRINT 180, (F MAG(K), K 
PRINT 170, (F ANG(K), K 
CONTINUE 


IF (K1 .EQ. 5) PUNCH 434, 
NR=NR-1 

TF (K2 (Sears GOmiOmS On 
IF (K2 .EQ. 4) K8 3 

IF (K2 .EQ. 3) K8 2 
Ki2tn=—0) 

KX = K2 + 1 

GO TO (561,551,551), KX 


* 57.29577951D0 


1,J) 


Ken 


180, 


(SRES(K), 


1,NN) 


K= 


1,NN) 


PRINT 533, NN,N, C(1), Z(2), C(2), 2(3), C(2), 2(4), C(4), 
SOURCE, RECVRS(40), FREQ 


FORMAT (1H1, 215, 10F10.4 
ICTR=0 

R LOS(1) 
LEVEL(1) 
DO 562 I 
P LEV(1)=40. 

J COU(1)=4 

J COUNT(1I)=-6 

CONTINUE 

IF((K2 .£Q. 2).AND.(NR .G 
GO TO 561 

NR = 5 

NL = NN 

PRINT 524, NL 

FORMAT (18, 13H MODES IN 
bee 4 

IF (K9 .GT. 0) NL = KQ 
READ 20, RANGE, FINAL R, 
IF (K8 .EQ. 3) PUNCH 30, 
IF (RANGE) 563,1,564 

NN = NN + 1 

READ 11, K1, K2, K3, K4, 
PRINT 11, K1, K2, K3, K4, 
GO TO 301 

IF (FINALR .LE. 0.) GO TO 
FINAL R = FINAL R + 1.D-3 
IF (RANGE .GE. FINALR) GO 


120. 


tomo 
Ss 


) 


T. 5))GO TO 772 


SUM) 


STEP R 


RANGE, FINAL R, STEP R 


K5, K6, 
KS, K6, 


550 


TO 561 


K7, 
K7, 


R ATTEN = RANGE * ATTEN - 9.9429946 
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K8, K9 


K8, 


K9 


521 
522 


523 


520 
536 


537 


534 


530 


535 


545 


903 
538 


IF (K1 .GT. 5) RATTEN = 0.D0O 

IF (ni wut. 5) GO TO 536 

IF (K7 .EQ. 2) RATTEN = 4.3429448D0 * DLOG( FREQ) 
IF (RANGE * DB LOSS(NL) .LT. 15.04) GO TO 522 
NE oe WEL = W 

DO 520 I = 1,NL 

IF (K7 .LT. 2) GO TO 523 

FMAG(I) = PHASE V(1I) 

GO TO 520 

TEMP RE = EIGEN(I) * RANGE 

TEMPIM = EIGENI(1) *RANGE 

CALL HZERO(TEMPRE, TEMPIM) 

HZERO2(I) = HZEROR 

IF (K7 .EQ. 0) GO TO 520 

FMAG(I) = HZEROR**2 + HZEROI**2 

HZER2I(1I) = HZEROI 


L=0 
DO 540 u = 1,NR 
FF = 0 
FFI = 0 


TEMP = 0.DO 

DO 530 I = 1,NL 

K=L+I1 

IF (Ki .LT. 6) GO TO 537 

TEMP = TEMP + UU(K) 

GO TO 530 

IF (K7 .EQ. 0) GO TO 534 

TEMP = TEMP + FMAG(I) * (UUI(K)**2 + UU(K) **2) 
GO TO 530 

TEMIM = UUI(K) * HZERO2(I) + UU(K) * HZER21(1) 
TEMRE = UU(K) * HZERO2(I) - UUI(K) * HZER21(1) 
FF = FF + TEMRE 

FFI = FFI + TEMIM 

CONTINUE 

IF (K1 .GT. 5) GO TO 535 

IF (K7 .GT. 0) GO TO 535 

TEMP = FF*¥*2 + FFI**2 

T RE = TEMP 

RX = -4.3429448 * ALOG(T RE) + R ATTEN 

R LOSS(J) = RX 

TF (K4 .LT. 2) GO TO 545 

T RE = -4.3429448 * ALOG(UU(K)**2 + UUI(K)**2) 
PRINT 170, RECVRS(JU), RX, T RE 

IF (K4 .NE. 3) GO TO 545 

CONTINUE 

tL = L + NN 

IF (K8 .LT. 2) GO TO 540 

IF (K8 .EQ. 3) GO TO 538 

LPCH = -RLOSS(J) * 10.D0 + 1400.5D0 

IF (LPCH .LT. 0) LPCH = 0 

IF (LPCH .GT. 999) LPCH = 999 

LOSPCH(JU,LL) = LPCH 

IRNG = RANGE / 1000.D0 


IF (LL .EQ. 25) PUNCH 903, IRNG, (LOSPCH(U,LLL),LLL 


FORMAT (15,2513) 
GO TO 540 
CONTINUE 
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1,25) 


684 
685 
686 
687 
688 
689 
690 
691 
692 
693 
694 
695 
696 
697 
698 
699 
700 
701 
702 
703 
704 
705 
706 
707 
708 
709 
710 
711 
712 
713 
714 
715 
716 
717 
718 
719 
720 
721 
722 
723 
724 
725 
726 
727 
728 
729 
730 
731 
732 
733 
734 
735 
736 
737 
738 
739 
740 


540 


712 


783 


776 


785 


787 


LOSS(J) 


= (1 
IF (LOSS(J) .LT. 
) eGo 


IF (LOSS(J 
CONTINUE 


40.05 - RX) * 10. 


0) LOSS(U) = 0 
999) LOSS(JU) = 999 


GO TO (770,780,716) ,KX 
C PLOT DB LOSSES 


COL(15)=1HI 
COL(39)=1HI 
COL(63)=1HI 
COL(87)=1HI 
COL(111)=1HI 
COL(27)=1HX 
COL(51)=1HX 
COL(75) =1HXx 
COL(99) =1HXx 
I PLACE = 135 


DO 787 J = 1,NR 
I PLACE = I PLACE - 24 
IPLOT = P LEV(J) 


IF (I PLOT .GT 


QoQ ¢A7 


af 


- R LOSS(J) 
0) GO TO 776 


P LEV(J) = P LEV(J) - 20. 

J COUNT(JU) = JU COUNT(J) - 2 

J COU(JU)=U COU(J)-2 

COL(I PLACE + 1) = 1HO 

IF (P LEV(J) - 100.) 778,779,781 


JC = J COU(JU) + 1 


COL(IPLACE) = ING(JC) 

GO TO 783 

COL(I PLACE) = 1HO 

GO TO 782 

JC = JU COUNT(J) + 1 
COL(IPLACE) = ING(UJC) 
COL(I PLACE - 1) = 1H1 
GO TO 783 

IF (I PLOT .LT. -9) GO TO 784 
GO TO 785 

P LEV(J) = P LEV(J) + 20. 


J COU(JU)=JU COU(J)+2 
JU COUNT(J) = JU COUNT(J) + 2 


GO TO 786 

IPP = I PLACE + IPLOT 
COL(IPP) = 1H+ 
CONTINUE 

GO TO 750 


C CONTOUR LOSS FIELD 
DO 590 JJ = 1,120 


780 
590 


620 


600 


COL(JJ) = 1H 
COL(31)=1HI 
COL(61)=1HI 
COL(91)=1HI 


DO 640 JJ = 2,41 


LEV = 1 

IF (RLOS (JJ) 
GO TO 610 
LEV=LEV+1 

GO TO 620 


eLT. 


CONTR(LEV)) GO TO 600 
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610 IF (LEV .EQ. LEVEL(JUJ)) GO TO 640 
650 IF (LEV .GT. LEVEL(JUJ)) GO TO 660 
GO TO 670 
660 II = LEVEL(JuJ) 
GO TO 680 
670 II = LEV 
680 JuJu = 124 - 3*Ju 
COL(JJJ) = J SMBL(IT) 
LEVEL(JJ) = LEV 
640 CONTINUE 
COL(1) = 1HI 
PRINT 261, (COL(1I1), I1 = 1,119) 
716 DO 690 JJ = 1,120 
690 COL(UJJ) = 1H 
ICTR = ICTR + 1 
IF (ICTR .EQ. 10) GO TO 700 


GO TO 714 
700 TEMP = (RANGE + 1.) / 10000. 
IND = TEMP 


COL(2) = ING(IND+1) 
TEMP1 = IND 
TEMP = f(TFMD — TEMP1) * 10. 
IND = TEMP 
COL(3) = ING( IND+1) 
TEMP1 = IND 
IND = (TEMP - TEMP1) * 10. 
COL(5) = ING(IND+1) 
COL(4)=1H. 
COL(6)=1HK 
COL(7)=1HY 
COL(8)=1HD 
COL(9) = 1HS 
ICTR=0 

714 GO TO (710,712) ,K2 

710 COL(31)=1HI 

COL(61)=1HI 
COL(91)=1HI 
DO 720 JJ = 1,40 
TEMP = LEVEL(Ju) 


TEMPI = O. 

830 IF (LEVEL(JJ) .GT. LEVEL(JUJ+1)) GO TO 730 
GO TO 740 

730 II = LEVEL(JJ) - 1 
KK = 1 


860 EX = (CONTR(II) - R LOS (JJ) )/ (RLOS (JUJU+1) - CONTR(ITI)) 
0O 760 LL = 1,3 
IF (EX .LT. TEST(LL)) GO TO 800 
760 CONTINUE 
LL = 4 
800 vJuJLL = 125 - 3*JU = LL 
COL(JUJLL) = JU SMBL(ITI) 
GO TO (810,820) ,KK 
810 LEVEL(JJ) = LEVEL(JJ) = 1 
GO TO 830 
740 IF (LEVEL(JJ) .LT. LEVEL(JUJ+1)) GO TO 840 
GO TO 720 
840 II = LEVEL(Ju) 
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798 
799 
800 
801 
802 
803 
804 
805 
806 
807 
808 
809 
B10 
81i 
812 
813 
814 
815 
816 
817 
818 


Cc 


KK=2 
GU IU BOU 
820 LEVEL(JJ) = LEVEL(JJ) + 1 
GO TO 740 
720 LEVEL(JJ) = TEMP 
COL(1) = 1HI 
750 PRINT 261, (COL(11), I1 = 1,119) 
261 FORMAT (1X, 119A1) 
GO TO 581 
PRINT DB LOSSES 
770 PRINT 580, RANGE, (R LOSS(K), K = 1,NR) 
Lue tb 24 
iP (Ub seu. 23) LL Se 4 
580 FORMAT (F9.0, 2X, 18F6.1) 
581 RANGE = RANGE + STEP R 
IF (K8 .NE. 3) GO TO 560 
PUNCH 980, (LOSS(I), I= 1,NR) 
980 FORMAT (2613) 
GO TO 560 
999 STOP 
END 
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25 


26 


30 


35 


SUBROUTINE SETUP 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION LAMBDA, LAMBDI 
COMMON /HAN/ H2R,H21,H1R,H11,H2PR,H2PI,H1PR,H1PI,EXPONT 
COMMON /EXPO/ EXSUM, CNTR, RATIO(25) 
COMMON/DETMNT/ A(25,4),Q(25,4) 
COMMON/INPUT/ Z(12), N, OMEGA, V, VI, CON(12), GSQ(12), 
1 CAY(12), LAMBDA, LAMBDI, G(12) 
2,RHO(12), GI(12), G SQI(12), CAYI(12) 

COMMON /LIMIT/ TLIM, EXPON , SLIM 
COMMON/PARTS/ ZT(12),ZTI(12),2B(12),ZBI(12) 
DENOM = v * V + VI * VI 

LAMBDA = OMEGA * V / DENOM 

LAMBDI = -OMEGA * VI/ DENOM 

M=N-1 

DO 10 I = 0,M 

IF (I .EQ. 0) GO TO 35 

TFEOZRENG— 74) GONT|ON 25 

IF (G(I) .LT. 0.) GO TO 25 

We lu) 8 (2M NS Zi) 2 22 

If (Ze (lh. =7.5) ZE = 1-75 

GO TO 26 

CONTINUE 

ZE = Gil) * (z(l+t) = Z(1)) + ZR 

IF (ZE .LT. SLIM) ZE = SLIM 

CONTINUE 

ZQ = GI(I) * (Z(I+1) - Z(1)) + ZI 
ZB(1) = ZE 

ZBI(1) = zQ 

CALL HANKEL(ZE,ZQ,0) 

ZB(1) = ZE 

ZBI(I) = ZQ 

RATIO(2+1) = EXPONT 

A(2*1,1) = H2R * RHO(I) 


Q(2*1,1) = H2I * RHO(1) 

A(2*1I,2) = H1R * RHO(I) 

Q(2*1I,2) = HII * RHO(1) 

A(2*I+1,1) = H2PR * G(I) - H2PI * GI(I) 
Q(2*I+1,1) = H2PI * G(I) + H2PR * GI(1) 
A(2*1I+1,2) = HIPR * G(I) - H1PI * GI(T) 
Q(2*1+1,2) = HIPI * G(I) + HIPR * GI(I) 
CONTINUE 


GSABS = G SQ(I+1)**2 + G SQI(I+1)**2 

X1 = CAY(I+1) - LAMBDA 

X2 = CAY(I+1) + LAMBDA 

X3 =CAYI(1I+1) - LAMBDI 

X4 =CAYI(I+1) + LAMBDI 

X = X1 * X2 -— X23 * X4 

Y= X1 * X4 + X3 * X2 

ZT(I+1) = (X * G SQ(I+1) + Y * G SQI(I+1)) / GSABS 
ZTI(1I+1) = (Y * G SQ(I+1) - X * G SQI(I+1)) / GSABS 
ZR = ZT(1+1) 

ZI = ZTI(I+1) 


ZE = ZR 
ZQ = ZI 


59 


36 
41 


40 


45 


10 


FaaGZ Rie Gili 


I 
S = cAy(i 
T = CAYI( 


+1) 
I+1) 


1,20 


-7.5) GO TO 40 


TEMP = S**2 + T**2 
TEMPI = (LAMBDI * S —- LAMBDA * T) / TEMP 


TEMP = (LAMBDA * S + LAMBDI * T) / TEMP 

ZR = ((1. + TEMP) * ( 1. - TEMP) + TEMPI**2) * CON(I+1) 
R= ZRey/ 75) 

IF (DABS(R-1.) .LT. 1.0-3) GO TO 41 

S = LAMBDA + (S - LAMBDA) / R 

ZI = -2. * TEMPI * TEMP * CON(I+1) 

74am) AR 

Ran EMD) Pal 

CONTINUE 

CALL HANKEL(ZR,ZI,0) 

ZT(1+1) = ZR 

7TT(T+1) = 7J 

RATIO(2+1+1) = EXPONT 

IF (I .NE. 0) GO TO 45 

A(1,3) = H2R * RHO(1) 

Q(1,3) = H2I1 * RHO(1) 

A(1,4) = H1R * RHO(1) 

Q(1,4) = H1I * RHO(1) 

GO TO 10 

CONTINUE 

A(2*I,3) =-H2R * RHO(I+1) 

Q(2*1,3) =-H2I * RHO(I+1) 

A(2*1,4) =-H1iR * RHO(I+1) 

Q(2*1,4) =-H1I * RHO(I+1) 

A(2*I+1,3) =-H2PR * G(I+1) + H2PI * GI(I+1) 
Q(2*1+1,3) =-H2PI * G(I+1) - H2PR * GI(I+1) 
A(2+#1+1,4) =-H1PR * G(I+1) + H1PI * GI(I+1) 
Q(2*1I+1,4) =-H1PI * G(I+1) - H1PR * GI(I+1) 
CONTINUE 

A( 2*N-2 ,4) = 0. 

Q( 2*N-2 ,4) = 0. 

A( 2*N-1 ,4) = 0 

RETURN 

END 
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(G(I1+1) * S + GI(I+1) * T) / (S**2 + T**2) 
/ CON**2 


ODNOUSWNH — 


30 


40 


50 


60 


70 


500 
80 


600 


SUBROUI INE VETNT(N,DET,DETI) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION DET, DETI 

COMMON /EXPO/ EXSUM, CNTR, RATIO(25) 
COMMON/DETMNT/ A(25,4), Q(25,4) 
DLOSS = 1. 

CNTR = 0. 

DET = A(1,3) 

DETI = Q(1,3) 

LIM =N+N-393 

DO 100 I[=1,LIM,2 


J=ztl 
K = 4 
su 
M=3 
1olet= a 
GO TO 500 
U= 1 +4 41 
K = 2 
Lead 
M=1 
GO TO 600 


Map 22 


or 


ie) 
I 
2 
O TO 500 
3 
3 


or 


n 
' 
=~ 


-EQ. LIM) GO TO 70 


= A(L,M)*A(L,M) + Q(L,M)*Q(L,M) 


DWODNDHYHROKHSHROKSROSCOKCC 
H 
bebo 


cI 


GO TO (10,50), II 
TD = A(uU,K) - (A(L,M)*B - Q(L,M)*BI) 
TOI = Q(uU,K) — (A(L,M)*BI + Q(L,M)*B) 
= TD*#2 + TDI**2 
TEMP = A(JU,K)**2 + Q(U,K)**2 
TEMP = TEM / TEMP 
IF (II .£Q. 2) GO TO 92 
IF (II .EQ. 4) GO TO 92 

= Q(U,K) * 10.D-18 
A(J,K) = A(J,K) * 10.D-18 
IF (TEMP .GT. 10.0-35) GO TO 92 
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(A(JU,K)*A(L,M) + Q(U,K)*Q(L,M)) / C 
= (Q(U,K)*A(L,M) -— A(U,K)*Q(L,M))/C 


100 


GO TO 90 

A(J,K) = TD 

Q(U,K) = TOI 

GO TO (700,40,60,70),II 

Cc = DET*A(U,K) - DETI*Q(U,K) 
DETI = DET*Q(JU,K) + DETI*A(U,K) 
DET =C 

GO TO (30,100), II 

CONTINUE 

RETURN 

END 
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OONOuUDSLWHD — 


102 


103 


10 


104 


105 


SUBROUTINE RCOEF (K3) 

IMPLICIT DOUBLE PRECISION (A-H,O-Z) 

COMMON/INPUT/ Z(12), N, OMEGA, Vv, VI, GCU(12), GSQ(12), 
1 CAY(12), LAMBOA, LAMBDI, G(12) 
2,RHO(12), G1(12), G SQI(12), CAYI(12) 

DIMENSION RR(12), RI(12), RA(12), RT(12), CYSQ(12), CYSQI(12) 
COMMON/REFL/ AF(12,200), AG(12,200), BF(12,200), BG(12,200), 
2 EIGEN(350), EIGENI(350),BR(25,4), BI(25,4), CB(12), CBI(12), 

3 CAYSQ(12), CAYSQI(12), NN 


NM = N - 1 

I = k3 

IF (I .GT. NM) I = NM 

vw S i a {i 

Ki =su eel 

IF (NN .NE. 1) GO TO 102 
Ls ta 4 


TEMP = CB(I)**2 + CBI(1)**2 
CY = OMEGA * CB(I) / TEMP 


CYI = -OMEGA * CBI(I) / TEMP 

CYSQ(I) = Cy**2 - CYI**2 

CYSOUKK) = CY © Cir 

CYSQI(I) = CYSQI(1) + CYSQI(1) 

EL SQ = C YSQ(I) - EIGEN(NN)**2 + EIGENI (NN) **2 
ELSQI = C YSQI(I) - 2.D0 * EIGEN(NN) * EIGENI(NN) 


TEMP = ELSQ + DSQRT (ELSQ*+*2 + ELSQI**2) 
IF (TEMP .LE. 0.D0) GO TO 107 

EL = DSQRT (TEMP * .5D0) 

ELI = ELSQI / (EL + EL) 

A = AF(I,NN)*BR(U,2) - AG(I,NN)*BI(uU,2) 
* + BF(I,NN)*BR(J,1) -— BG(I,NN)*BI(UJ,1) 

B = AF(I,NN)*BI(JU,2) + AG(1,NN)*BR(J, 2) 
* + BF(I,NN)*BI(U,1) + BG(1I,NN)*BR(U,1) 

E = AF(I,NN)*BR(K,2) - AG(I,NN)*BI(K,2) 
* + BF(1,NN)*BR(K,1) — BG(1I,NN)*BI(K,1) 

F = AF(I,NN)*BI(K,2) + AG(1,NN)*BR(K,2) 
* + BF(I,NN)*BI(K,1) + BG(1,NN)*BR(K,1) 

C = (F * EL - E * ELI) / (ELSQ + ELSQI) 
D=-(E * EL + F * ELI) / (ELSQ + ELSQI) 
TEMP = (A + C)**2 + (B + D)**2 
RR(I) = (A**2 -— Cx*2 + Be*2 - D**2) / TEMP 


RI(1) -2.D0 * (A * D- B * C) / TEMP 
FORMAT (10013.5) 

RA(I) = 0 

IF (CB(I) .GT. V) GO TO 104 

RX = CB(I) / V 

RA(I) = ACOS(RX) * 57.296 

RT(1I) = RR(I)**2 + RI(1)**2 

RT(I) = 1.D0 / RT(I) 

RI(I) = -DATAN2 (RI(I), RR(I)) * 57.296D0 
RR(I) = -4.34294D0 * DLOG (RT(I)) 

IF (K3 .NE.1) GO TO 108 

Tare ly at 

IF (I .LT. N) GO TO 110 

CONTINUE 

PRINT 106, (RR(I), I = 1,NM) 

PRINT 106, (RI(I), I = 1,NM) 
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106 


108 
109 


107 


PRINT 106, (RT(I), I 
PRINT 106, (RA(I), I 
FORMAT (9G13.4) 
RETURN 


PRINT 109, Z(L), RR(I), RT(1), RI(I), 


FORMAT (9H AT DEPTH, 


1,NM) 
1,NM) 


F7.0,8H YD, 


R 


RA(1) 
7F9.4 


»7H DB, OR,D12.4, 


* 8H, PH A = ,F9.3,16H DEGREES, GR A =,F8.2,9H DEGREES. ) 


RETURN 
EL = 0.D0 


ELI = DSQRT ( DABS (ELSQ) ) 


GO TO 103 
END 
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ODNOuUDBWH— 


1 


GFONMONOWDFOODNOUSWNH— * 


GTOANMONDPFPOMDNIOUHBWNDA — * 


SUBROUTINE HANKEL (ZR, ZI, IH) 
COMMON /HAN/ H2R,H21I,H1R,H11,H2PR,H2PI,H1PR,H1IPI,R 
INTEGER FLPS, FLQUAD 

DOUBLE PRECISION ZMLA2, TLIM, R 
COMMON /LIMIT/ TLIM,EXPONT 


DOUBLE PRECISION 
A,B0,B,CO,C,D0O,D, c4, 


ZR,Z1,H2R,H21,H1R,H1L1,H2PR,H2PI,HIPR,HIPI,AO, 
D4,K,K1,K2,CON4,STORE1,STORE2, 


STORE3,STORE4, STORES, STOREG6, STORE7, STORES, STOREO, STOR10,STORI1, 
STOR12,STOR13,STOR14,STOR1I5,STOR16,STOR1@,C1,C2,C3,CPR,CPI,CTHR, 


CTP,FR,FI,FPR,FPI, 


FIR,F1L1,F2R,F21,GR,GI,GPR,GPI,GIiR, 


G11,G2R,G21,H11R,H111,H12R,H121,H11PR,H11P1,H12PR,H12PI, 


H21PR,H21PI,H22PR,H22PI ,H21R,H211,H22R,H221, 


SPI,S2, 


EXPONT ,ABK, 


ZRT41, 


XPR(40), 


THR,X,XR,XI,XPR,XPI,YR,YI, 


PI,SR,SI,SPR, 
ZM, ZMSQ, ZM1R,ZM11, 


ZRTM2R,ZRIM21,ZRTM4R,ZRIM41,ZRT2R,ZRT2ZI,ZRT2M,ZRI4GR, 


Z32R,Z2321,STP,STOR17,STHR,C5,D05,FP112 
DIMENSION A(40),8(40),C(40) ,D(40) ,C4(20) ,D4(20),ZMLA2(40), 


DATA (A(I),1=1,36) / 


-1.5507278615487157D-001, 


XPI1(40), C5(20), D5(20), ZMLA5(20) 


5.16909287182905237D-03, 


-7.17929565531812830D-05, 
-2.58993349758951237D0-09, 
-2.01519879986734545D-14, 
-5.20045934975470047D-20, 
-5.66054875234532879D-26, 
-3.03137585006604588D-32, 
-8.89081245106713441D-39, 
-1.54547567290139313D-45, 
-1.69172458673478018D-52, 
-1.22347235365465573D-59, 
-6.07825428023425146D0-67, 
-2.142372753117290340-74, 
-5.50471327313964415D-82, 
-1.05526034966989126D-89, 
-1.53977168218007537D-97, 
-1.74020422875830830D-105, 
-1.54687790339646370D-113, 


DATA (B(1),1=1,36) / 
-5.6524893762022989D-002, 


-1.49536755984187802D-05, 
-3.99403728590245199D-10, 
-2.52780770480649350D-15, 
-5.57276830865629078D-21, 
-5.34066309073303316D-27, 
-2.570196682611954820-33, 
-6.87508809232765398D-40, 
-1.10221782500302268D-46, 
-1.12255630004727929D-53, 
-7.60687925589328600D-61, 
-3.56156318721341813D-68, 
-1.18880450319548524D-75, 
-2.90462369773880309D-83, 
-5.31361078500669380D-91, 
-7.421557145296762240-99, 
-8.05038914195299455D-107, 
-6.88468878345390325D-115, 


5 .438860344937975960-07, 
8. 46383495944285089D-12, 
3.65072246352779973D-17, 
5.97753948247666720D-23, 
4.49249900979787999D-29, 
1.76038086531129261D-35, 
3.94096296589855249D-42, 
5.39998488085741835D-49, 
4.77888301337508527D-56, 
2.85191690828591078D-63, 
1.18901687798009614D-70, 
3.56705420099448941D-78, 
7.89545793623012643D-86, 
1.317428651273272490-93, 
1.68834614274131072D-101, 
1.68919067050893836D-109, 
1.33859285513712678D-117/ 


9.58568948616588476N-0RA. 
1.16784715962060000D-12, 
4.21301284134415583D-18, 
5.99222398780246320D-24, 
4.00950682487464952D-30, 
1.42314323511182437D-36, 
2.92308167190801615D-43, 
3.71117112795630532D-50, 
3.06709371597617291D-57, 
1.72023501942408096D-64, 
6.77618566821426585D-72, 
1.92925106003811301D-79, 
4.06810041700112477D-87, 
6.48792525641842955D-95, 
7.96988525053346460D-103, 
7.66265861598419431D-111, 
5.84835948 305632284D0-119/ 
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1.34583080385769022D-03, 


TO NMVNDPOMDNINODWHSWH — * 


GTONMVGNWDFODNOUDBWND— * 


HUH HHH HH KH 


eH HK 


DATA ( 


-6. 
Alo 
-8. 
0 
A} 5 
CU 
ales 
-2. 
26 
—\q 
-8. 
=2. 
Ol 
S'\9 
malite 
ithe 
"Wo 


DATA ( 


6 
=O) 
—3)5 
Oi} 
Si) 6 
Pio 
-3. 
=O}. 
-6. 
-4. 
=2. 
Ole 
-2. 
-4. 
=6). 
-8. 
C5 


DATA(C4(1), 


-. 1826444261 146441383D003/ 
DATA(D4(I), 


C(1),1=1,36) 


-3.1014557230974314D-002, 


52663241392557118D-06, 
52349029269971316D-10, 
76173391246671935D-16, 
79326184474300016D-21, 
61729964352723680D-27, 
3935996 34307428970-34, 
891662223631305190-40, 
91599183566300591D-47, 
86732980802505116D-54, 
88226515946870112D0-61, 
56092152145669220D0-69, 
78230227677570174D-76, 
632184666433306210-84, 
18568578614594524D-91, 
62081229703165829D-99, 
72297448391911713D-107, 
44568028354809692D-115, 


MET =A", 26) if 


-2.2609957504809195D-001, 


49536755984187802D-04, 
39045965744392318D-09, 
56117895057428569D-14, 
56037512642376142D-19, 
81582545084923127D0-25, 
02807867304478193D-31, 
16254052247072083D-38, 
73153269001571794D-45, 
51082654027421986D-52, 
86840272377170304D-59, 
49309423104939269D-66, 
03491422428568780D-74, 
38179143214581853D-81, 
67597749080589054D-89, 
97626371657895650D-97, 
05038914195299455D-105, 
29777011046113744D-113, 


Wel oI) 7 

.1041666666666666666D000, 
-2290716053934337712D001, 
-9062847663874030839D001, 
.2032967817611733257D002, 
.3609376712592949187D002, 
-5635611849738394099D002, 
-8111724483308463849D002, 
-1103774469890957246D003, 
- 1441391353710093869D003, 


Read lO) 7 
. 1458333333333333333D000, 


-.2190010740010626122D001, 
-.8910269876251375731D001, 
-.2013810597032928847D002, 
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3.88490024638426856D-08, 
4.23191747972142544D-13, 
1.40412402443376913D-18, 
1.86798108827395850D-24, 
1.18223658152575789D-30, 
4.00086550298021048D-37, 
7.881925931797104970-44, 
9.64283014438824705D-51, 
7.70787582802433108D-58, 
4.19399545336153351D-65, 
1.60677956483796775D-72, 
4.45881775124311176D-80, 
9.180765042128053990-88, 
1.43198766442747010D-95, 
1.72280218647072522D-103, 
1.62422179856628689D-111, 
1.21690259557920616D-119/ 


1.24613963320156502D-06, 
2.21890960327913999D-11, 
1.05325321033603896D-16, 
1.85758943621876359D-22, 
1.48351752520362032D-28, 
6.11951591098084481D-35, 
1.43231001923492791D-41, 
2.04114412037596793D-48, 
1.87092716674546548D-55, 
1.15255746301413424D-62, 
4.94661553779641407D-70, 
1.52410833743010928D-77, 
3.45788535445095606D-85, 
5.90401198334077089D-93, 
7.73078869301746066D-101, 
7.89253837446372014D-109, 
6.37471183653139190D-117/ 


6.46136608978631547D-04, 


9.42081562700383155D-03, 


- .5876374421296296294D000, 


-5115246914604383039D0001, 
-1413420435039637896D002, 
-2764948541118776109D002, 
-4566262114618547916D002, 
-6817431262327333036D002, 
-9518494701182174927D0002, 
- 1266948822584689090D003, 
- 1627327073751970376D003, 


-5242693865740740734D000, 
- 4986 228080782250417D001, 


-.1396107513192384956D002, 
-.2744104880120119211D002, 


-.3586970305443466737D002, 
-.5610363044805757000D002, 
-.8083919802992951438D002, 
-.1100759893558574910D003, 
—.1438158255226323270D003, 
-.1826430810500566861D003/ 
DATA(C5(1), 1=1,19) / 
-.8020833333333333322D000, 
-.3778635359389885240D001, 
-.6771926170404923857D001, 
-.9768433620097970692D001, 
-.1276620743800279721D002, 
-. 15764630887 78327385D002, 
-.1876343942052581175D002, 
-.2176249668603070326D002, 
-.2476104285001339985D002, 
-.2772097924164600240D002/ 
DATA(D5(1), I=1,19) / 
-.6770833333333333322D000, 
-.3712590308861503947D001, 
-.6721592615807936173D001, 
—.9726176955678393129D001, 
-.1272907520499096141D002, 
-.1573111965670700007D002, 
-.1873266140557309258D002, 
-. 21733874940 4898291 0D002, 
—.2473404635295374440D002, 
—.2758686258346552664D002/ 


ee He HH HK 


ee RH HHH HH HH 


H * He HH HR HH KR 


DATA (ZMLAS(1I), I=4#,17) / 1.E9,715., 
* 22.6, 18.5, 16.6, 14.7, 14., 12.9, 


DATA LA2, LAS /36,17/ 
DATA (ZMLA2(1),1=1,40) / 


-.4542393171194671533D002, 
- .6790874395729851569D002, 
= .9489495581082038481D002, 
-.1263823305425793323D003, 
-.1624125356051020867D003, 


-.2285545023696682453D001, 
-.5274623160711059862D001, 
- .82699549413402740750001, 
-.1126721374768006806D002, 
-.1426535892246475774D002, 
-.1726399730395179650D002, 
- . 20262943441 40785724D002, 
— .2326204842605378820D002, 
- .2625430055269059493D002, 


- 2202914798206278015D001, 
-5218010289498877822D001, 
-8224186533352991837D001, 
-11227766967105/87810U02, 
-1423017624892113601D002, 
-17231939820272207020002, 
- 2023330232819351955D002, 
- 2323435992394561375D002, 
- 2622252877088035571D002, 


207., 103., 47., 36.4, 27., 


W265 ViseS5 Vaso Goe / 


12.6944301D-12,4.7348244D-6,7.08037130-4,9.6398932D-3, 
24.9271494D-2,1.5267301D-1,3.5324772D-1 ,6.7835277D-1, 
31.1475215D0,1.7731141D0,2.5617749D0,3.9957181D0,5.2159327D0, 
46.6020702D0,8.1490972D0,9.8514112D0,1.2320405D1,1.4917923D1, 
51.7163634D1,1.9540309D1 ,2.303024601,2.6446844D1 ,2.929282001, 


63.3549744D1,3. 
75.3917166D1,5. 
88.0163399D1,8. 
OY .OS7/S1 Sebe, I 
DATA Cl £ Wo 
DATA C2 7 Oo 


754124601 ,4. 
9582285D1,6. 
6945766D1,9. 
1039828D2,1. 


0802980D1,4 
4537858D1,7 
2594255D1,9 
1820169D2/ 


57735 026918962576D0 / 
66666 666666666666D0 / 


C3 / 0.86602 540378443864D0 / 

PI / 3.14159265358979324D0 / 
FPI12 / 1.30899693899574718D0/ 

CON4 / .7071067811865475244D0/ 
9.304367 16929229427D-01 


BO = 6.78298725144275871D-01 
CO = 4.65218358464614714D-01 
DO = 6.78298725144275871D-01 

K = 0.85366721883895156D0 
ZMSQ = ZR*ZR +ZI1*ZI 

RZR = ZR 

TEMP = ZI 

IF (TEMP .LT. O.) TEMP = -TEMP 
IF (RZR .LE. 0.) GO TO 51 
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-5784933D1 ,5.0287133D1, 
-0701377D1,7.5978472D1, 
-983446D1, 


171 IF (RZR .GT. 4.4) GO TO 120 


172 TEM1 = 7. — .2632 * RZR«*2 
173 IF (TEMP .GT. TEM1) GO TO 120 
174 GO TO 53 

175 51 IF (RZR .LT. -9.) GO TO 120 
176 TEM1 = 4.4 + .1375 * RZR 

177 IF (TEMP .GT. TEM1) GO TO 120 
178 53  FLPS = 1 

179 STORE1 = ZR*ZR-ZI¥*ZI 

180 STORE2 = 2.*ZR*ZI 

181 XR = STORE1*ZR -STORE2*ZI 

182 XI = STORE1#ZI +STORE2*ZR 

183 DO 55 MLS=1,LA2 

184 IF (ZMSQ - ZMLA2(MLS)) 62,62,55 
185 55 CONTINUE : 

186 62 FR = AO 

187 FI = 0.0 

188 XPR(1) = XR 

189 XPI(1) = XI 

190 DO 65 M = 1,MLS 

191 FR=FQ+A(M)*XPR(M) 

192 FI=FI+A(M)*XPI(M) 

193 XPR(M+1)=XR*XPR(M)—-XI*XPI(M) 
194 XPI(M+1)=XI*XPR(M)+XR*XPI (M) 
195 65 CONTINUE 

196 GR=BO 

197 GI=0.0 

198 DO 72 M = 1,MLS 

199 GR=GR+B(M)+*XPR(M) 

200 GI=GI+B(M)+*XPI(M) 

201 72 CONTINUE 

202 X =ZR+GR-ZI*GI 

203 GI=ZR*GI+Z1*GR 

204 GR=X 

205 SR=-C1*(GI-FI-FI) 

206 S1=C1+(GR-FR-FR) 

207 H2R=GR-SR 

208 H21=GI-SI 

209 GO TO 317 

210 120 FLPS = 0 

211 ZM = DSQRT(ZMSQ) 

212 ZRT2M = DSQRT(ZM) 

213 IF (ZR .LT. 0.D0) GO TO 125 
214 ZRT2R = DSQRT (0.5D0 * (ZR + ZM)) 
215 ZRT2I = ZI / (ZRT2R + ZRT2R) 
216 Z32R = ZR*ZRT2R - ZI*ZRT21 
217 Z321 = ZR¥ZRT2I + ZI*ZRT2R 
218 GO TO 130 

219 125 ZRT2I = DSQRT (0.5D0 * (ZM - ZR)) 
220 IF (ZI .LT. 0.D0) ZRT2I = -ZRT2I 
221 ZRT2R = ZI / (ZRT2I + ZRT21) 
222 Z32R = ZR*ZRT2R - ZI*ZRT21 
223 Z321 = ZR*ZRT2I + ZI*ZRT2R 
224 ZM1R = DABS(Z321) 

225 IF (ZM1R .LT. TLIM) GO TO 130 
226 R = (TLIM / ZMIR) 

227 Z32R = Z32R * R 
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228 R = DCBRT(R) 


229 ZRT2R = ZRI2R * R 

230 ZRT2I = ZRT2I * R 

231 ZRT2M = ZRT2M * R 

232 R=R* R 

233 ZM = ZM * R 

234 ZMSQ = ZM**2 

235 ZR = ZR * R 

236 ZI = ZY * R 

237 130 ZRT4R = DSQRT (0.5D0 * (ZRT2R + ZRT2M)) 
238 ZRT41I = 0.500 * ZRT2I / ZRT4R 
239 ZRTIM4R = ZRT4R/ZRT2M 

240 ZRTM4I = -ZRT41/ZRT2M 

241 IF (ZR .GT. 0.) GO TO 210 

242 IF (ZMiR .LT. TLIM) GO TO 210 
243 ABK = ABS(K2) 

244 IF (2321 .GT. 0.) GO TO 205 
245 K1 = K * EXPONT 

246 K2 = K / EXPONT 

247 Z32I = -TLIM 

248 GO TO 220 

249 905 «2 = K * FXPONT 

250 K1 = K / EXPONT 

251 Z32I = TLIM 

252 GO TO 220 

253 210 K2 = C2 * Z32I 

254 S2 = DEXP(K2) 

255 K2 = K*S2 

256 K1 = K/S2 

257 220 THR = FPI12 - C2 * Z32R 

258 STHR =DSIN(THR) 

259 CTHR =DCOS(THR) 

260 STP = -C3*CTHR +0.5*STHR 

261 CTP = -C3*STHR -0.5*CTHR 

262 TEMP = DABS (Z32R) 

263 TEM1 = DABS (2321) 

264 IF (TEMP .LT. TEM1) TEMP = TEM1 
265 230 DO 235 ML = 1,LA5 

266 IF (TEMP .GT. ZMLAS(ML)) GO TO 250 
267 235 CONTINUE 

268 250 CONTINUE 

269 VR = 2321 

270 YI = -Z32R 

271 CALL CFR (YR, YI, F2R, F2I, C4, C5, ML) 
272 CPR = F2R 

273 CPI = F2I 

274 STORE3=K2* (ZRTM4R* F2R-ZRIM41*FQI ) 
275 STORE4=K2*(ZRTIM41+F2R+ZRIM4R*F QI ) 
276 H22R =STORE3*CTRR-STORE4*STHR 
277 H221 =STORE3*STHR+STORE4*CTHR 
278 IF (ZR) 280,270,270 

279 270 FLQUAD =0 

280 GO TO 300 

281 280 IF (ZI) 290,310,310 

282 290 FLQUAD = 1 

283 300 H2R = H22R 

284 H2I = H22I1 
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350 


92 


94 


380 


GO TO 317 

FLQUAL = -) 

YR = -Z321 

YI = Z32R 

CALL CFR (YR, YI, FIR, FII, C4, C5, ML) 
CPR = FIR 

CPI = F1I 


STORE5=K1*(ZRTIM4R*F1R-ZRIM41*F11) 
STORE6=K1¥*(ZRTIM4R*F1I+ZRTM41*F1R) 
H21R=STORES*CTP-STOREG*STP 
H211=STORES*STP+STORE6*CTP 
H2R=H21R+H22R 

H21=H211+H221 

IF (IH .EQ. 2)GO TO 80 

IF (FLPS .NE. 1) GO TO 320 

H1R = GR+SR 

H11I = GI+SI 

GO TO 362 

IF (FLQUAD .LT. 0)GO TO 340 

YR = -Z321 

YI = Z32R 


CNG CER (WR Yito FUR Fil, C45 @Sq Liud)) 


STORE7=K1*(ZRTIM4R*F1R-ZRIM41*F11) 
STORE8=K1*(ZRTIM41*F1R+ZRIM4R¥*F11) 
H11R=STORE7*CTHR+STORE8*STHR 
H111=STORE7*(-STHR)+STORE8*CTHR 
IF (FLQUAD .LE. 0) GO TO 360 
STORE9=K2*(ZRTM4R*F2AR-ZRIM4GI*FQ2I) 
STOR10=K2*(ZRTIM41*F2R+ZRIM4R¥FQ21 ) 
H12R = STOREQ*CTP+STOR10*STP 

H12I = STORE9*(-STP)+STOR10*CTP 


H1R = H11R+H12R 
H1I = H111+H121 
GO TO 362 
H1R = H11R 
HiIl = H111 


IF (IH .EQ. 1)GO TO 999 

IF (FLPS .NE. 1) GO TO 380 
FPR = CO 

FPI = 0.0 

DO 92 M = 1,MLS 

FPR=FPR+C(M) *XPR(M) 
FPI=FPI+C(M)*XPI(M) 

X =-(STORE1*FPR-STORE2*FPI) 
FPI=-(STORE1*FP1I+STORE2*FPR) 


FPR = X 
GPR = DO 
GPI = 0.0 


DO 94 M = 1,MLS 
GPR=GPR+D(M) *XPR(M) 
GPI=GP1+D(M)*XPI(M) 
SPR=-C1+*(GPI-FPI-FPI) 
SPI=C1*(GPR-FPR-FPR) 
H2PR=GPR-SPR 
H2PI=GPI-SPI 

GO TO 414 

YR = 2321 
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@PRT,S 


J.CFR 


390 
400 


410 


999 


YI = -Z32R 


CALL CFR (YR, YI, G2R, G21, D4, DS, ML) 


STOR11 = K2*(ZRT4R*G2R-ZRT41*G21 ) 
STOR12 = K2*(ZRT4R*G2I1+ZRT41*G2R) 
H22PR=STOR11*STHR+STOR12*CTHR 
H22PI=STOR11*(-CTHR) +STOR12*STHR 
IF (FLQUAD .LT. 0) GO TO 410 


H2PR = H22PR 
H2PI = H22PI 
GO TO 414 
YR = -Z32I1 
YI = Z32R 


CALL CFR (YR, YI, GIR, G1I, D4, DS, ML) 


STOR13 = K1*(ZRT4R*G1R-ZRT41*G11) 
STOR14 = K1*(ZRT4R*G11+ZRT41*G1R) 
H21PR=STOR13*(-STP) -STOR14*CTP 
H21PI=STOR14*(-STP) +STOR13*CTP 
H2PR = H21PR+H22PR 

H2PI = H21PI+H22PI1 

IF (IH .EQ. 2) GO TO 999 

ie (SLRS olileo 1) CO) ue) Geo 

H1PR = GPR+SPR 

HiPI = GPI+SPI 

GO TO 999 

IF (FLQUAD .LT. 0) GO TO 440 

YR = -Z32I1 

YI = Z32R 


CALL CFR (YR, YI, GIR, Gil, D4, D5, ML) 
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STOR15 = K1*(ZRT4R*GIR -ZRT41*G11) 
STOR16 = K1¥*(ZRT4R*GII +ZRT41*G1R) 
H11PR = STOR15*STHR -STOR16*CTHR 
H11PI = STOR15*CTHR +STOR16*STHR 
IF (FLQUAD .GT. 0) GO TO 470 

HiPR = H11PR 

HiPI = H11PI 

GO TO 999 

STOR17 = K2*(ZRT4R*G2R -ZRT41*G21) 
STOR18 = K2*(ZRT4R*G2I +ZRT41*G2R) 
H12PR = STOR17*(-STP) +STOR18*CTP 
H12PI = STOR17*(-CTP) -STOR18*STP 
H12PR+H11PR 

H1PI = H12PI+H11PI 

CONTINUE 

RETURN 

END 


=x 
= 
vu 
a 
ot 
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ODNOUSBAWA— 


SUBROUFINE CFR(X, Y, SR, SI, A, B, M) 
IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
DIMENSION A(1), B(1) 

SR = 0.D0 

SI = 0.D0 

DO 10 J = 1,M 

I=M-J5+41 

TEMR X + SR + B(1) 

TEMI Yc Si 

TEMP = A(I) / (TEMR**2 + TEMI**2) 
SR = TEMR * TEMP 


SI = -TEMI * TEMP 
CONTINUE 

SR =-SR + 1.00 
RETURN 

END 
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10 
11 


SUBROUTINE HZERO(PARTR, PARTI) 
IMPLICIT DOUBLE PRECISION(A-H,O-Z) 
DOUBLE PRECISION PARTI, PARTR, K2, IMZ12 
COMMON /AHZERO/ HZEROR,HZEROI 

D2 = PARTR**2 + PARTI**2 

K2 .7978845608+*EXP( PARTI ) 

D = SQRT(D2) 

RLZ12 = (SQRT((PARTR + D)/2.))/D 

IF(D - PARTR) 9,9,10 

IMZ12 = 0 

GO TO 11 

IMZ12 = (SQRT((D - PARTR)/2.))/D 
COST = COS(PARTR) 

SINT =-SIN(PARTR) 

HZEROR = K2*(RLZ12*COST - IMZ12*SINT) 
HZEROI = K2*(IMZ12*COST + RLZ12*SINT) 
RETURN 

END 
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APPENDIX B: SAMPLE RUN 


This appendix gives a brief discussion of the input-output, then lists an input deck 
and shows the resulting output. The input deck is really three separate runs that are stacked 
to run consecutively. The input to a single run consists of several parts given in table B1. 
The table gives the number of cards and the location of the FORTRAN input statements in 
Program MAIN. The last three of these are open-ended. That is, more modes, receiver 
depths, or ranges are read in until a blank card specifies the end of the list. A blank range 
card sends the program back to the beginning. A negative range sends the program back to 
read a new source and new receivers after reading another key card. The program halts 
when a blank “‘n and freq” card is encountered. 


Table B2 gives most of the functions of the key card by which integers are read into 
control keys 1-9. Some of these will require additional information, which is read in imme- 
diately following the key card. 


The output of the program is usually printed through FORTRAN print statements. 
Cards are also punched when key 5 is 10 or key 8 is 2, 3, or 4. In the first case each card 
contains a complete eigenvalue that can be read into future runs. 


When key 8 = 2, propagation losses for 25 consecutive ranges per card are punched 
for each receiver depth, with a maximum of 5 receiver depths. The losses can be read into a 
plot program with a format of (5X,25F3.1). Each loss must then be subtracted from 140. 
This format allows losses to tenths of a dB from 40.1 to 140.0 dB. 


When key 8 is 3, losses for up to 26 receiver depths are punched on one card for 
each range. These cards can be used in a contour plotting program. They can be read witha 
format of (26F3.1) and must also be subtracted from 140. 


Table B1. Input cards to the normal mode program. 


Location 
in Program 
Input Function Number of Cards MAIN 
Control keys Selects options 1 or more 37-65 
nand freq Determines number of layers and fre- 1 66 
quency; also halts program 
Profile Specifies depths, sound speeds, gradients, 71-85 
attenuations, and densities 
Modes Searches for or specifies modes 1 or more 224 
Source and receivers | Specifies a source depth and one or 2 or more 461-463 
more receiver depths 
Ranges Specifies a sequence of ranges; also 1 or more 616 


directs continuation 
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Table B2. Functions of the control keys. 


Setting Effect Function Affected 


No output Depth functions 


Print 
Plot on printer 
2 0 Print losses Propagation losses 
1 Contour on printer 
3 0 No output Reflection coefficients 
1 Print all interfaces 
k>1 Print interface k 
4 k>0 Change levels and symbols Contour on printer 
5 0 Sum only those given Number of modes 
1 Add to existing sum 
k+10 Punch modes on cards 
6 0 Long print Steps in mode iteration 
1 Short print 
2 Shortest print 
i 0 Phased addition Mode sum 
1 Random-phase addition 
8 1 Change T-lim 
Punch losses for up to 5 Loss plot input 
receivers 
3 Punch losses for up to 26 Contour plot input 
receivers 
9 0 No effect Number of modes 
k Use only Ist k modes 


The first profile in the input-output is a surface duct, 100 m deep. For the 500 Hz 
frequency, 3 modes are found by searching from a phase velocity of 1520.5-1523 m/s. Two 
additional modes are found by extrapolation. Forty receiver depths are then specified from 
3 to 120 m, and propagation loss contours are drawn for a source at a depth of 20 m. The 
modes are added in random phase, and loss contours of 80, 90, and 100 dB are requested to 
be represented by the symbols 8, 9, and 0. A negative range then causes the program to read 
new control keys, source and receivers. The depth functions are then printed out as ampli- 
tudes and phase angles and propagation losses are computed. 


The second profile consists of two negative gradients over two layers of sediments in 
shallow water. A velocity discontinuity exists at the top of each sediment layer. Negative 
numbers in the input for the attenuation at the bottom of the sediment layers serve as flags 
to request that the gradients at the top of the layers have no imaginary parts and that the 
attenuation at the bottom will be whatever results from this. The change in ImC from 37.9 
to 23.7 in the deeper sediment layer indicates that the attenuation changed by about 60 
percent through the layer. 


Wd 


A final layer of negative sound speed gradient must always be added. A gradient of 
-0.1 is chosen here for the top of this layer. 


The first three modes are determined by reading in approximate values. The magni- 
tudes of the depth functions are plotted on a log scale at 2-m depths from 30 to 80 m. 
Reflection coefficients are computed at interface 2, which is the water-sediment interface. 


The final profile is a deep-water profile with a 40-m deep sediment layer. The atten- 
uation increases from 2 dB/km to 2.5 dB/km through this layer. The first mode, the first 
bottom-reflected mode, and a higher bottom-reflected mode are found. Each step of the 
mode iteration is printed out. Reflection coefficients are again computed. The amplitudes 
and phase of the depth functions are printed out at 500-m depth and at each even 1000-m 
depth for a 100-m source depth. 


On the last two profiles, a much larger set of modes is required to compute correct 
propagation losses. 


The sample run given here required 6 seconds on a UNIVAC 1110, Exec 8 operating 
system. The cost of the run was $1.20. 
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NORMAL*MODE(0) . INPUT 


INPUT DECK STARTS AT LINE 3, ENDS AT LINE 65. 


' 
2 123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 
3 1 1 Tate ace 
4 100. 90. 80. -1000. 
5 0 
6 098 
7 2 500. 
8 100. 
9 1520. 
10 0 
1 .017 -.1 
12 0 
13 0 
14 0 
15 -1520.5 1523. 30 
16 -1. 2 
17 0 
18 10. 
19 3. 120. ay 
20 0 
21 4000. 100000. 4000. 
22 -1. 
23 1 
24 10. 
25 30. 120. 30 
26 0 
27 5000. 100000. 5000. 
28 0 
29 2 2 1 
30 5 1500. 
31 O. 51. TB. 1308) 373.3 
32 1542.2 1536.8 1606.45 1684 
33 1523.42 
34 1.5 1.5 -.1 
35 a2 .73 73 
36 =1. -1. 
37 1.68 1.91 1.91 
38 1527.18 16 
39 1530.64 13 
40 1533.49 11 
a1 0 
42 60. 
43 30. 80. 2. 
44 0 
45 ry) 
46 1 6 
47 8 100. 
48 55. 146. 402. 960. 2286. 4390. 4430. 
49 1544.9 1542.6 1517.9 1495.0 1483.2 1497.8 1541.7 
50 1533.4 
51 J -.1 
52 02 .025 
53 .025 
54 1.54 2.5 
55 -1483.5 1484.5 10 
56 -1533.4 a 1534.4 10 


Uy 


-1600. 2 1602. 10 


1000. 5000. 1000. 


123456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789 
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APPENDIX C: HANKEL FUNCTION PARAMETERS 


This appendix gives the FORTRAN statements for two programs associated with the 
modified Hankel functions. Program PWRTRN computes the power series coefficients, dy), 
from eq (57), then determines the truncation points from eq (59). The truncation points 
for the other three sets of coefficients can be determined by changing line 9. Different com- 
puter word lengths can be accommodated by changing line 16. 


The second program, CFC, determines the asymptotic series coefficients C,, from 


eq (72), then determines the continued fraction coefficients as indicated by eq (81)-(83). 
The second set of coefficients can be determined by changing the 4 in line 11 toa 16. 


89 


ODNOW AWAD — 


C ** THIS PROGRAM DETERMINS TRUNCATION POINTS FOR THE POWER SERIES. 


50 


60 


30 


20 


40 
10 


PROGRAM PWRTRN 


IMPLICIT DOUBLE PRECISION(A-H,O0-Z) 


DIMENSION D(50), 


mCi) 20s 
ALOGD(1) = 0. 
P = 3. 


ALOGD(59) 


ALOGD(I) = ALOG10(D(I)) 


CONTINUE 
PRINT 60, D, 


ALOGD 


FORMAT (10E12.6) 


DH = 18. 


K 
- K 
ALOGD(K ) 
LOGD(M) 


1 


QZrHPHNvV 
(o) 
" 
= 
ROSZPPU-Z20- 


wort 
fo) 
= 


M- 1 


AZ 


- ALOGD(M) + DH) / 3. / P 
-GT. -1.1) GO TO 20 


- ALOGD(M+1) - 3. 


AZSQ = AZ * AZ 


PRINT 40, L, 
FORMAT (215, 
CONTINUE 
END 


MM, Z, 
4E15.8) 


-GT. 0.) GO TO 20 


EXP (Z * 2.3025851) 


AZ, AZSQ 


* Z 
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PROGRAM CFC 

C ** THIS PROGRAM COMPUTES A SET OF SERIES COEFFICIENTS AND THEN 

C ** COMPUTES THE CORRESPONDING CONTINUED FRACTION COEFFICIENTS. 
IMPLICIT DOUBLE PRECISION(A-H,O-Z) 
DIMENSION COEF(21,23,3), CHECK(20), C(82), S(10), A(20), B(20) 


ODNOUAWH = 


20 


100 


110 


120 


150 


130 


6) S Vo 

BOTTOM = 1 

TOP = 1 

po 21 = 1,45 

X = 48 * I 
Y=9* (I + I1- 1)**2 - 4 


CCM) 2 CU) oY 7 2 
CONTINUE 

PRINT 20, (C(I), I = 1,40) 
FORMAT (5G20.9) 

FORMAT (/) 

DO 100 I = 1,1 
COEF(I,1,3) = 
COEF(I,I+1,3) 
COEF(I,1+2,3) 
CONTINUE 

A(1) = C(2) 
COEF(2,2,3) = 1. 
DO 140 I = 3 
DO 110 uU = 2 
COEF(I,u,1) 
COEF(I,U,2) 
COEF(I,uU,3) 
CONTINUE 

IF (I .EQ. 3) GO TO 150 
CON = 0. 
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K=I1- 3 

DO 120 JU = 3,1 

K + 1 

= C(K) * COEF(I,U-1,3) + CON 
AT = C(K) * COEF(I,U-1,2) + AT 
BT = C(K) * COEF(I,U-1,1) + BT 
CONTINUE 

PRINT 160, CON, AT, BT 
CHECK(1-2) = BT 

A(I-2) = -(CON + C(K+1)) / AT 
CONTINUE 

CCN = O. 


K =K + 1 

CON = C(K) * COEF(I,U-1,3) + CON 
AT = C(K) * COEF(I,U-1,2) + AT 
BT = C(K) * COEF(I,U-1,1) + BT 
CONTINUE 

PRINT 160, CON, AT, BT 

PRINT 11 
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140 


160 


30 
200 


B(I-2) = -(C 
DO 140 JU = 2 


COEF(I,JU,3) = COEF(I,U,3) + A(I-2) * COEF(I,u,2) + B(I-2) * 


* COEF(I,UJ,1) 
CONTINUE 
PRINT 20, A, B, CHECK 
FORMAT (5G20.9) 


K = -2 

J=0 

DO 30 M = 1,18,3 
J=uJd+ 3 

Is & iS sp) 83 


PUNCH 200, (A(I), I 
PUNCH 200, (B(I), 1 
CONTINUE 


FORMAT (5X, 1H*, 3(£21.15, 


END 


ON + A(I-2) * AT + C(K+1)) / BT 
I 


)) 
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APPENDIX D: MODE FOLLOWER PROGRAM IN FORTRAN 
The FORTRAN statements of the Mode Follower program are given here. This is 


the main body of the program. The following auxiliary subroutines from appendix A are 
required: SETUP, DET, HANKEL, and CFR. 
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OODNOUAWNH— 


PROGRAM MFOLLO 
IMPLICIT DOUBLE PRECISION (A-H,0O-Z) 


COMMON /INPUT/ Z(10),N,OMEGA,V,VI,GCU(10),GSQ(10) ,CAY(10),LAMBDA,L 


1AMBDI,G(10),RHO(10),GI(10) ,GSQI(10),CAYI(10) 

COMMON /DETMNT/ A( 21,4) ,Q(21,4) 

REAL INCA, INCB, INCC, INCD, INCE, LAMBDA, LAMBDI 

DIMENSION T(4), PV(4), W(8), WI(8), CB(10), CBI(10), C(10), 


1 CAY SQ(10), GAMMA(10), DPK(10), GCUI(10), CI(10), CR(10),PVI(4) 


2 , CAYSQI(10), SR(4), S1(4) 
CHNG = 1. / 8192. 


CHNGI = 0. 
4 CONTINUE 
Crt nO - TOTAL STEF LIMIT, K1, K2 PRINT KEYS, K3 = i nEEPS SAME PRUriLeE 
C** FOR NEXT RUN. 


READ 10, KO, K1, K2, K6, K3, TLIM, BLIM, RATIO, EX 
PRINT 10,KO, K1, K2, K6, K3, TLIM, BLIM, RATIO, EX 
10 FORMAT (514, 4610.1) ' 
IF (TLIM .EG. 0.) TLIM = 
IF (BLIM .EQ. 0.) BLIM = 
IF (EX .EQ. 0.) EX = 28. 
RLIM = 10.**EX 
IF (RATIO .EQ. 0.) RATIO = 2. 
IF (KO .EQ. 0) KO = 300 
IF (K3 .NE. 0) GO TO 128 
30 READ 1240, N,FREQ  ,ATTEN 
C** STOP IF N = 0. THIS IS THE ONLY PROGRAMED STOP. 
IF (N.EQ.G) GO TO 1200 
PRINT 1250, N,FREQ 
C** PARAMETERS READ IN BELOW ARE THOSE AT THE TOP OF EACH LAYER. 
C** READ IN VELOCITIES. 
READ 1260, (C(I),I=1,N) 
PRINT 1280, (C(I),I=1,N) 
C** READ IN DEPTHS. 
READ 1260, (Z(I),I= 
PRINT 1280, (Z(I),I 
C** READ IN GRADIENTS 
READ 1260, (GAMMA(I),I=1,N) 
PRINT 1280, (GAMMA(I),I=1,N) 
C** READ IN ATTENUATION FACTOR IN LOSS PER KILOMETER. 
READ 1260, (DPK(I),I=1,N) 
PRINT 1280, (DPK(I),I=1,N) 
C** READ IN DENSITIES (BLANK INPUT IMPLIES SEA WATER DENSITY). 
READ 1260, (RHO(I),I=1,N) 
PRINT 1280, (RHO(I),I=1,N) 
128 CONTINUE 
NUMBER = 1 
UX = 0 
C** NX = VARIABLE, NY = LAYER NUMBER, NZ = CONTINUITY 
READ 119, NX, NY, NZ, PK, VALL,DP, V, VI, STEP, STEPI 
PRINT 21, NX, NY, NZ, PK, VALL,DP, V, VI, STEP, STEPI 
119 FORMAT (312, 4X, 7010.2) 
21 FORMAT (10H VARIABLE ,I2, 10H LAYER NO , 12,12H CONTINIITTV 
* 12, / 7G15.5) 
PK = PK — DP 
C** START NEW CYCLE BY INCREMENTING VARIABLE. 
107 PK = PK + DP 


1,N) 
=1,N) 
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IF (DP) 108,999,109 
C** CHECK IF DESIRED LIMIT OF VARIABLE HAS BEEN REACHED. 
108 IF (PK .LT. VALL) GO TO 3 
GO TO 133 
109 IF (PK .GT. VALL) GO TO 3 
133. GO TO (131,101,102,103,104,105) ,NX 
131 FREQ = PK 
GO TO 106 
101 C(NY) = PK 
IF (NZ .NE. 0) GO TO 106 
134 IF (NY .EQ. N) GO TO 135 
GAMMA(NY) = O. 
IF (NY .LT. 2) GO TO 106 
135 GAMMA(NY-1) = O. 
GO TO 106 
102 Z(NY) = PK 
IF (NZ .EQ. 1) GO TO 106 
IF (NY .LT. N) GO TO 134 
IF (NUMBER .EQ. 1) GO TO 106 
C(NY) = O. 
GO TO 106 
103 GAMMA(NY) = PK 
IF (NZ .NE.0) GO TO 106 
J = NY + 1 
DO 121 I. = uJ,N 
(i) S Wo 
121 CONTINUE 
104 DPK(NY) = PK 
GO TO 106 
105 RHO(NY) = PK 
106 CONTINUE 


C#* COmMPLETG PROFILE ** 
DO 100 I=1,N 
C** SET UNSPECIFIED DENSITIES TO 1.02 (SEA WATER). 
IF (RHO(I).NE.O.) GO TO 40 
RHO(1I)=1.02 
40 IF (1.EQ.1) GO TO 50 
C** COMPUTE VELOCITY AT BOTTOM OF PREVIOUS LAYER. 
TEMP=CI(I-1)**2 
TEMDR=C(I-1)**2 
TEMDI=(TEMDR+TEMDR+TEMDR-TEMP) *CI( 1-1) 
TEMDR=(TEMDR-TEMP-TEMP-TEMP) *C(I-1) 
TEMP=(GAMMA(I-1)+GAMMA(1-1))*(Z(1)-Z(1I-1) )-C(I-1) 
TEMDEN=TEMP**2+CI(I-1)**2 
TEM1=(TEMDI*CI(I-1)-TEMDR*TEMP) /TEMDEN 
TEM11=(-TEMDI*TEMP-TEMDR*CI(I-1))/TEMDEN 
CB(1)=SQRT(.5*(TEM1+SQRT(TEM1**2+TEM11**2) ) ) 
CBI(1I)=TEM11/(CB(I)+CB(1)) 


50 IF (C(1).NE.0) GO TO 60 

C** IF VELOCITY WAS UNSPECIFIED USE VELOCITY AT BOTTOM OF PREVIOUS LAYER 
C(1)=CB(1) 

60 IF (DPK(1).NE.0.) GO TO 70 
CI(1)=0. 
GO TO 80 


C** IF ATTENUATION IS TO BE APPLIED TO A LAYER, COMPUTE COMPLEX VELOCITY 
C** KEEP ABSOLUTE C EQUAL TO GIVEN REAL C FOR SIMPLICITY. 
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70 


80 


C** 


110 


Cre 


120 


C** 


71 


250 


C*x 


TEMP=545 74%. * FREQ 
TEMDI=DPK(I1)*C(1) 
TEMDR=TEMP**2+TEMDI**2 
CI(1)=TEMDI*TEMP*C(1I)/TEMDR 
C(1)=TEMP**2*C(I)/TEMDR 


IF (GAMMA(I).NE.O.) GO TU 100 


IF (1.EQ.N) GO TO 9O 
COMPUTE GRADIENT IF NOT GIVEN. 


GAMMA(1)=(C(I+1)**2-C(1)**2)*C(1)/(2.*C(14+1)**2*(Z(1+1)-Z(1))) 


i (US BOG) Kel) 8), Se) 
GO TO 100 


REDUCE LAYERS BY ONE IF FINAL POINT ONLY DEFINES GRADIENT IN LAST LAYER. 


N=N-1 
CONTINUE 


COMPUTE USEFULL QUANTIES ** 
OMEGA=6.283185307+FREQ 
DO 120 1=1,N 
TEMP=C(1I)**2+CI(1)**2 
CAY(1)=OMEGA*C(1)/TEMP 
CAYI(1)=-OMEGA+CI(1)/TEMP 


CAYSQ(1)=CAY(1)**2-CAYI(1)**2 


CAYSQI(1)=2.*CAY(1)*CAYI(1) 
TEMDR=-2.+*GAMMA(1)*CAYSQ(1) 
TEMDI=-2.*GAMMA(I) *CAYSQI(1) 


GCU(1I)=(TEMDR*C(1)+TEMDI*CI(1))/TEMP 
GCUI(1)=(TEMDI*C(1I)-TEMDR*CI(1))/TEMP 
TEMP=EXP(ALOG(GCU(1)**2+GCUI(1)**2)/6.) 


GI(1)=TEMP*SIN(ATAN2(GCUI(1),ABS(GCU(I)))/3.) 


G(1)=SQRT(TEMP**2-GI(1)**2) 


IF (GAMMA(I).LT.O.) GO TO 110 


G(1)=-G(T) 
GI(I)=-GI(1) 


vm TS A LAYER STRENGH PARAMETER USED ONLY TO 


XMI=~-GI(1)*(Z(1I+1)-Z(1) ) 
XM=-G(1)*(Z(1I+1)-Z(1)) 
GSQI(1)=2.*G(1)*GI(1) 
GSQ(1)=G(1)**2-G1(1)**2 
IF (JX .GT. 0) GQ TO 113 


felel 
eur 


Wine 


GO TO INITIAL 3 STEPS OR TO THE STANDARD STEP. 


IF (NUMBER - 4) 71,111,122 
CALL SETUP 
CALL DETNT(N,DET,DETI) 
VEL=DET 
VELI=DETI 
DELTA=STEP 
DELTI=STEPI 
IF(DELTA.NE.O.)GO TO 250 
IF(DELTI.EQ.0. )DELTA=.01 
SIZE2=100. 


Liven meen ween 
Wain WIIen wwe 


IF (K6.LT.3) PRINT 1320, V,VI,DET,DETI,A(21,4),Q(21,4) 


ITERATE FOR MODE UP TO 7 STEPS. 
DO 310 J=1,12 
V=V+DELTA 
VI=VI+DELTI 


DO NOT PERMIT IMAGINARY PART TO BECOME NEGATIVE. 


IF (VI) 260,270,280 
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260 
270 


C** SET UP DETERMINANT FOR PHASE VELOCITY V + VI 


280 


300 


125 


C** DISCONTINUE ITERATION AFTER 2ND STEP IF CORRECTION STEP IS MORE THAN 


DELTI=DELTI-VI 
VI=1.E-18 


CALL SETUP 
CALL DETNT (N,DET,DETI) 
IF (K6.NE.1) GO TO 300 
PRINT 1330, V, VI, DET, DETI, 
TEMNR = DELTA 
TEMNI = DELTI 
TEMDR=VEL-DET 
TEMDI=VELI-DETI 


TEMDEN=TEMDR* TEMDR+T EMDI*TEMDI 


IF (TEMDEN.EQ.0.) GO TO 32 


Toimainu-7TomNR*TEMDR+TEMNI*TEMDI 
TEMINU=TEMNI*TEMDR-TEMNR*TEMDI 


SLR = TEMRNU / TEMDEN 

SLI = TEMINU / TEMDEN 

IF (J .GT. 3) GO TO 125 

SR(4-NUMBER) = SLR 

SI(4-NUMBER) = SLI 

DEWAN= DETee*s SER DET I = SET 

DELTI = DET * SLI + DETI * SLR 
SIZE=DELTA*DELTA+DELTI*DEL 


C** PREVIOUS STEP. 


310 


320 


51 


111 


122 


112 
113 


IF ((SIZE.GT.SIZE2).AND.(JU.GT.2)) GO TO 320 


SIZE2=SIZE*2. 
VEL=DET 
VELI=DETI 
CONT INUE 
CONTINUE 
PV (4-NUMBER) Vv 
PVI(4-NUMBER)= VI 
NUMBER = NUMBER + 1 


GO TO 107 
C*x* START STANDARD STEP, 
INCA = DP 
INCB = DP 
INCC = DP 
INCD = -INCB —- INCC 
T(1) = INCB * INCD 
T(2) = INCB * INCC 
T(3) = INCD * INCC 
DO 112 IS = 1,3 
W(IS+4) = -SR(IS) / T(IS) 
WI(IS + 4) = -SI(IS) / T(IS) 
W(IS) = -PV(IS) / T(IS) 
WI(IS) = -PVI(IS) / T(IS) 
INCD = INCA + INCB 
INCE = INCD + INCC 
T(1) = INCD * INCE 
T(2) = INCA * INCE 
T(3) = INCA * INCD 
SLOP = 0. 
SLOPI = 0. 
SUM = 0. 
SUMI = 0. 


SLR, 


0 


TI 
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EXTRAPOLATE PHASE VELOCITY AND SLOPE. 


00 114 IS = 1,3 

SLOP = SLOP + W(IS + 4) * T(IS) 
SLOPI = SLOPI + WI(IS+4) * T(IS) 
SUM = SUM + W(IS) * T(IS) 

SUMI = SUMI + WI(IS) * T(1IS) 


Vv = SUM 
VI = SUMI 
CALL SETUP 


CALL DETNT (N,DET,DETI) 


C** EVALUATE DETERMINANT AT THE EXTRAPOLATED POINT. 


VEL = DET 
VELI = DETI 


C*x* ITERATE FOR THE ROOT USING EXTRAPOLATED SLOPE. 


DET * SLOP = DETI * SLOPI 

DET * SLOPI + DETI * SLOP 

IF (K1 .EQ. 1) PRINT 1330, V, VI, DET, DETI, DELTA, 
V = V + DELTA 

VI = VI + DELTI 

IF (VI .GE. 0.) GO TO 124 
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DEER T= aD pales Vil 
CHNGI = CHNGI - VI 
VI = 0 


C** RE-EVALUATE AT NEW POINT. 


124 


C** 


123 


CALL SETUP 
CALL DETNT (N, DET, DETI) 
TEMNR = DELTA 
TEMNI = DELTI 
TEMDR=VEL-DET 
TEMDI=VELI-DETI 
TEMDEN=TEMDR* TEMDR+TEMDI*TEMDI 
IF (TEMDEN .EQ. 0.) GO TO 123 
TEMRNU=TEMNR*TEMDR+TEMNI*TEMDI 
icM I NU=7EMNI*TEMDR-TEMNR*TEMDI 


EVALUATE SLOPE (RECIPROCAL ACTUALLY USED). 


SLR = TEMRNU / TEMDEN 

SLI = TEMINU / TEMDEN 

DELTA Bey 62 SUR — Piste 2 Sibi 

DELTI DET * SLI + DETI * SLR 

IF (K1 .£Q. 1) PRINT 1330, V, VI, DET, DETI, DELTA, 


CORRECT PHASE VELOCITY TO BEST VALUE. 


Vie= Ve > DELTA 
Wil a Wat ap pyelE rae 
TEMP = Vx*2 / (TEMNR**2 + TEMNI**2) 


IF (TEMP .LT. RLIM) GO TO 123 
IF (TEMP .LT. 1.€34) GO TC 141 
SLR = SLOP 


IF NOT, USE EXTRAPOLATED SLOPE. 


SLI = SLOPI 
GO TO 141 
CONTINUE 


C**x IF SO, FIND 1 - RATIO OF SLOPES. 


TEMDEN = (SIR**2 + SLI¥*2) 

TEMDR = SLR * SLOP + SLI * SLOPI - TEMDEN 
TEMDI = SLR * SLOPI - SLI * SLOP 

TEMP = (TEMDR**2 + TEMDI*¥*2) / TEMDEN**2 
IF (TEMP .GT. TLIM) GO TO 116 


98 


DELTI 


DELTI 


WAS INCREMENT LARGE ENOUGH TO PERMIT EVALUATION OF SLOPE. 


C** SLOPE RATIO TOO GOOD. DOUBLE STEP. 


141 DP = DP * RATIO 
GO TO 117 


116 IF (TEMP .LT. BLIM) GO TO 117 


PRINT 130, PK,V,VI,DET,DETI,SLR,SLI,TEMP,OBLOSS ,NUMBER 
130 FORMAT (1X,E14.6,£16.9,£13.7,6E10.3,15) 
C** SLOPE RATIO TOO POOR. HALVE STEP. 


IF (NUMBER .LT. 7) GO TO 126 
PK = PK — DP 

DP = DP / RATIO 

INCA = DP 

UX = JX + 1 


IF (K2 .EQ. 1) PRINT 118, PK, 


C** STOP ON 5 SUCCESSIVE FAILURES. 
IF (UX .LT. 5) GO TO 107 
PRINT 810, N, FREQ 
810 FORMAT (14, G12.5) 
3 DO 801 I = 1,N 


PRINT 800, C(I), Z(I), GAMMA(I), DPK(I), RHO(I), 


800 FORMAT (10G12.5) 
801 CONTINUE 

GO TO 4 
126 PRINT 127 , N, TEMP 


127 FORMAT (7H NUMBER, I3, 22H FAILED, SLOPE RATIO 


V, VI, DET, DETI 
MODE IS LOST. 


C** UPDATE ALL QUANTITIES FOR NEXT STEP. 


117 INCC = INCB 
INCB = INCA 
INCA = DP 


PV(3) = PV(2) 
PVI(3) = PVI(2) 
Pv(2) = PV(1) 
PVI(2) = PVI(1) 


Pv(1) = V 

PVI(i) - vi 

UX = 0 

DENOM = v * V + VI * VI 
LAMBDI = -OMEGA * VI / DENOM 
DB LOSS = -8686. * LAMBDI 
SR(3) = SR(2) 

SR(2) = SR(1) 

SR(1) = SLR 

SI(3) = SI(2) 

SI(2) = SI(1) 

SI(1) = SLI 


GV = Vx¥*2 / (V - FREQ * (V-Pv(2)) / INCB) 


PRINT 118, PK,V,VI,DET,DETI,SLR,SLI, TEMP,OBLOSS,GV,NUMBER 


118 FORMAT (£15.7,£16.9,€13.7,6E10.3,F11.5,15) 


NUMBER = NUMBER + 1 
C** CHECK TOTAL NUMBER OF STEPS. 
IF (NUMBER .GT. KO) GO TO 3 
GO TO 107 
999 STOP 


1250 FORMAT (13,8H LAYERS,,F10.1,3H HZ) 


1260 FORMAT(SE10.4) 


1270 FORMAT (8H ATTEN =,G10.5,5HDB/KM) 


1280 FORMAT (8F14.5) 
1320 FORMAT (/,6E18.9) 
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G(I) 


,F10.6) 


342 
343 
344 
345 
346 
347 


1330 FORMAT (6E18.9) 


1240 FORMAT(1I2,E10.1, £10.2) 


1290 FORMAT (7X,6H RE M 


,8X,6H IM M ,8X,6H L/KYD,8X,6H RE C ,8X,6H IMC 


1 ,5X,12H RE C BOTTOM,4X,12H IM C BOTTOM) 


1200 CONTINUE 
END 
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